####################################################################
## author:    Lukas P. Fesenfeld
## contact:   lukas.fesenfeld@ir.gess.ethz.ch
## file name: fesenfeld et al_models_final.R
## Context:   The role and limits of strategic framing for promoting sustainable consumption and policy. 
## Date:      2019-11-02
## Summary:   Statistical models for framing experiment
##############################################

rm(list = ls())

# Set Working directory ---------------------------------------------------

# Load Packages ---------------------------------------------------
library("ggplot2")
library ("vctrs")
library ("ggstance")
library ("RColorBrewer")
library("sparsereg")
library("estimatr")


#Load Data--------

df <- read.csv("C:/Users/lf20q095/Downloads/data_clean_final_2019-11-28.csv", sep = ",", header = T)


#Datasubsets--------


df$country <- factor(df$country, levels = c("China", "Germany", "USA"))
df_comp <- subset(df, treatment_check_food =="True" | treatment_check_mobility =="True")

df_food_USA <- subset(df, df$country=="USA" & df$group =="Food")
df_food_Germany <- subset(df, df$country=="Germany" & df$group =="Food")
df_food_China <- subset(df, df$country=="China" & df$group =="Food")

df_food <- subset(df, df$group =="Food")

df_mobility <- subset(df, df$group =="Mobility")

df_mobility_USA <- subset(df, df$country=="USA" & df$group =="Mobility")
df_mobility_Germany <- subset(df, df$country=="Germany" & df$group =="Mobility")
df_mobility_China <- subset(df, df$country=="China" & df$group =="Mobility")

pd <- position_dodge(0.1)



#Balance tests------

m.gender_group <- glm(gender=="female" ~ group, df, family=binomial())
summary(m.gender_group)

m.gender_food_usa <- glm(gender=="female" ~ Food_, data= subset (df, country =="USA" &group =="Food" ), family=binomial())
summary(m.gender_food_usa)

m.gender_mobility_usa <- glm(gender=="female" ~ Transport_, data= subset (df,country =="USA" & group =="Mobility" ), family=binomial())
summary(m.gender_mobility_usa)

m.gender_food_germany <- glm(gender=="female" ~ Food_, data= subset (df, country =="Germany" &group =="Food" ), family=binomial())
summary(m.gender_food_germany)

m.gender_mobility_germany <- glm(gender=="female" ~ Transport_, data= subset (df,country =="Germany" & group =="Mobility" ), family=binomial())
summary(m.gender_mobility_germany)

m.gender_food_china <- glm(gender=="female" ~ Food_, data= subset (df, country =="China" &group =="Food" ), family=binomial())
summary(m.gender_food_china)

m.gender_mobility_china <- glm(gender=="female" ~ Transport_, data= subset (df,country =="China" & group =="Mobility" ), family=binomial())
summary(m.gender_mobility_china)

library(stargazer)
BTgender <- stargazer(m.gender_food_usa, m.gender_mobility_usa, m.gender_food_germany, m.gender_mobility_germany, m.gender_food_china, m.gender_mobility_china,
                                           type="html", 
                                           single.row = TRUE,
                                           title="Balance Test Gender", 
                                           align=TRUE,
                                           star.char = c("+", "*", "**", "***"),
                                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                                           digits = 3,
                                           column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                                           object.names = FALSE,
                                           notes.append = FALSE, notes.align = "l",
                                           notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                                           out='BTgender.html')

m.age_group <- lm(age ~ group, df)
summary(m.age_group)

m.age_food_usa <- lm(age ~ Food_, data= subset (df, country =="USA" &group =="Food" ))
summary(m.age_food_usa)

m.age_mobility_usa <- lm(age ~ Transport_, data= subset (df,country =="USA" & group =="Mobility" ))
summary(m.age_mobility_usa)

m.age_food_germany <- lm(age ~ Food_, data= subset (df, country =="Germany" &group =="Food" ))
summary(m.age_food_germany)

m.age_mobility_germany <- lm(age ~ Transport_, data= subset (df,country =="Germany" & group =="Mobility" ))
summary(m.age_mobility_germany)

m.age_food_china <- lm(age ~ Food_, data= subset (df, country =="China" & group =="Food" ))
summary(m.age_food_china)

m.age_mobility_china <- lm(age ~ Transport_, data= subset (df,country =="China" & group =="Mobility" ))
summary(m.age_mobility_china)

BTage <- stargazer(m.age_food_usa, m.age_mobility_usa, m.age_food_germany, m.age_mobility_germany, m.age_food_china, m.age_mobility_china,
                      type="html", 
                      single.row = TRUE,
                      title="Balance Test Age", 
                      align=TRUE,
                      star.char = c("+", "*", "**", "***"),
                      star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                      digits = 3,
                      column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                      object.names = FALSE,
                      notes.append = FALSE, notes.align = "l",
                      notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                      out='BTage.html')

m.educLow_group <- glm(educationGRP=="Low" ~ group, df, family=binomial())
summary(m.educLow_group)

m.educLow_food_usa <- glm(educationGRP=="Low" ~ Food_, data= subset (df,country =="USA" & group =="Food"), family=binomial())
summary(m.educLow_food_usa)

m.educLow_mobility_usa <- glm(educationGRP=="Low" ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ), family=binomial())
summary(m.educLow_mobility_usa)

m.educLow_food_germany <- glm(educationGRP=="Low" ~ Food_, data= subset (df,country =="Germany" & group =="Food"), family=binomial())
summary(m.educLow_food_germany)

m.educLow_mobility_germany <- glm(educationGRP=="Low" ~ Transport_, data= subset (df, country =="Germany" &group =="Mobility" ), family=binomial())
summary(m.educLow_mobility_germany)

m.educLow_food_china <- glm(educationGRP=="Low" ~ Food_, data= subset (df,country =="China" & group =="Food"), family=binomial())
summary(m.educLow_food_china)

m.educLow_mobility_china <- glm(educationGRP=="Low" ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ), family=binomial())
summary(m.educLow_mobility_china)

BTeducLow <- stargazer(m.educLow_food_usa, m.educLow_mobility_usa, m.educLow_food_germany, m.educLow_mobility_germany, m.educLow_food_china, m.educLow_mobility_china,
                   type="html", 
                   single.row = TRUE,
                   title="Balance Test Education Low", 
                   align=TRUE,
                   star.char = c("+", "*", "**", "***"),
                   star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                   digits = 3,
                   column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                   object.names = FALSE,
                   notes.append = FALSE, notes.align = "l",
                   notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                   out='BTeducLow.html')

m.educMedium_group <- glm(educationGRP=="Medium" ~ group, df, family=binomial())
summary(m.educMedium_group)

m.educMedium_food_usa <- glm(educationGRP=="Medium" ~ Food_, data= subset (df, country =="USA" & group =="Food" ), family=binomial())
summary(m.educMedium_food_usa)

m.educMedium_mobility_usa <- glm(educationGRP=="Medium" ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ), family=binomial())
summary(m.educMedium_mobility_usa)

m.educMedium_food_germany <- glm(educationGRP=="Medium" ~ Food_, data= subset (df, country =="Germany" & group =="Food" ), family=binomial())
summary(m.educMedium_food_germany)

m.educMedium_mobility_germany <- glm(educationGRP=="Medium" ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility" ), family=binomial())
summary(m.educMedium_mobility_germany)

m.educMedium_food_china <- glm(educationGRP=="Medium" ~ Food_, data= subset (df, country =="China" & group =="Food" ), family=binomial())
summary(m.educMedium_food_china)

m.educMedium_mobility_china <- glm(educationGRP=="Medium" ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ), family=binomial())
summary(m.educMedium_mobility_china)

BTeducMedium <- stargazer(m.educMedium_food_usa, m.educMedium_mobility_usa, m.educMedium_food_germany, m.educMedium_mobility_germany, m.educMedium_food_china, m.educMedium_mobility_china,
                       type="html", 
                       single.row = TRUE,
                       title="Balance Test Education Medium", 
                       align=TRUE,
                       star.char = c("+", "*", "**", "***"),
                       star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                       digits = 3,
                       column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                       object.names = FALSE,
                       notes.append = FALSE, notes.align = "l",
                       notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                       out='BTeducMedium.html')

m.educHigh_group <- glm(educationGRP=="High" ~ group, df, family=binomial())
summary(m.educHigh_group)

m.educHigh_food_usa <- glm(educationGRP=="High" ~ Food_, data= subset (df, country =="USA" & group =="Food" ), family=binomial())
summary(m.educHigh_food_usa)

m.educHigh_mobility_usa <- glm(educationGRP=="High" ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ), family=binomial())
summary(m.educHigh_mobility_usa)

m.educHigh_food_germany <- glm(educationGRP=="High" ~ Food_, data= subset (df, country =="Germany" & group =="Food" ), family=binomial())
summary(m.educHigh_food_germany)

m.educHigh_mobility_germany <- glm(educationGRP=="High" ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility" ), family=binomial())
summary(m.educHigh_mobility_germany)

m.educHigh_food_china <- glm(educationGRP=="High" ~ Food_, data= subset (df, country =="China" & group =="Food" ), family=binomial())
summary(m.educHigh_food_china)

m.educHigh_mobility_china <- glm(educationGRP=="High" ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ), family=binomial())
summary(m.educHigh_mobility_china)

BTeducHigh <- stargazer(m.educHigh_food_usa, m.educHigh_mobility_usa, m.educHigh_food_germany, m.educHigh_mobility_germany, m.educHigh_food_china, m.educHigh_mobility_china,
                          type="html", 
                          single.row = TRUE,
                          title="Balance Test Education High", 
                          align=TRUE,
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          digits = 3,
                          column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                          object.names = FALSE,
                          notes.append = FALSE, notes.align = "l",
                          notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                          out='BTeducHigh.html')

m.inc_group <- lm(inc_n ~ group, df)
summary(m.inc_group)

m.inc_food_usa <- lm(inc_n ~ Food_, data= subset (df, country =="USA" & group =="Food" ))
summary(m.inc_food_usa)

m.inc_mobility_usa <- lm(inc_n ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ))
summary(m.inc_mobility_usa)

m.inc_food_germany <- lm(inc_n ~ Food_, data= subset (df, country =="Germany" & group =="Food" ))
summary(m.inc_food_germany)

m.inc_mobility_germany <- lm(inc_n ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility" ))
summary(m.inc_mobility_germany)

m.inc_food_china <- lm(inc_n ~ Food_, data= subset (df, country =="China" & group =="Food" ))
summary(m.inc_food_china)

m.inc_mobility_china <- lm(inc_n ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ))
summary(m.inc_mobility_china)

BTinc <- stargazer(m.inc_food_usa, m.inc_mobility_usa, m.inc_food_germany, m.inc_mobility_germany, m.inc_food_china, m.inc_mobility_china,
                        type="html", 
                        single.row = TRUE,
                        title="Balance Test Income", 
                        align=TRUE,
                        star.char = c("+", "*", "**", "***"),
                        star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                        digits = 3,
                        column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                        object.names = FALSE,
                        notes.append = FALSE, notes.align = "l",
                        notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                        out='BTinc.html')

m.occupationWorking_group <- glm(occupationGRP=="Working" ~ group, df, family=binomial())
summary(m.occupationWorking_group)

m.occupationWorking_food_usa <- glm(occupationGRP=="Working" ~ Food_, data= subset (df,country =="USA" & group =="Food" ), family=binomial())
summary(m.occupationWorking_food_usa)

m.occupationWorking_mobility_usa <- glm(occupationGRP=="Working" ~ Transport_, data= subset (df,country =="USA" & group =="Mobility" ), family=binomial())
summary(m.occupationWorking_mobility_usa)

m.occupationWorking_food_germany <- glm(occupationGRP=="Working" ~ Food_, data= subset (df,country =="Germany" & group =="Food" ), family=binomial())
summary(m.occupationWorking_food_germany)

m.occupationWorking_mobility_germany <- glm(occupationGRP=="Working" ~ Transport_, data= subset (df,country =="Germany" & group =="Mobility" ), family=binomial())
summary(m.occupationWorking_mobility_germany)

m.occupationWorking_food_china <- glm(occupationGRP=="Working" ~ Food_, data= subset (df,country =="China" & group =="Food" ), family=binomial())
summary(m.occupationWorking_food_china)

m.occupationWorking_mobility_china <- glm(occupationGRP=="Working" ~ Transport_, data= subset (df,country =="China" & group =="Mobility" ), family=binomial())
summary(m.occupationWorking_mobility_china)

BToccupationWorking <- stargazer(m.occupationWorking_food_usa, m.occupationWorking_mobility_usa, m.occupationWorking_food_germany, m.occupationWorking_mobility_germany, m.occupationWorking_food_china, m.occupationWorking_mobility_china,
                   type="html", 
                   single.row = TRUE,
                   title="Balance Test Working", 
                   align=TRUE,
                   star.char = c("+", "*", "**", "***"),
                   star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                   digits = 3,
                   column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                   object.names = FALSE,
                   notes.append = FALSE, notes.align = "l",
                   notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                   out='BToccupationWorking.html')

m.occupationNotWorking_group <- glm(occupationGRP=="No paid working" ~ group, df, family=binomial())
summary(m.occupationNotWorking_group)

m.occupationNotWorking_food_usa <- glm(occupationGRP=="No paid working" ~ Food_, data= subset (df, country =="USA" & group =="Food"), family=binomial())
summary(m.occupationNotWorking_food_usa)

m.occupationNotWorking_mobility_usa <- glm(occupationGRP=="No paid working" ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ), family=binomial())
summary(m.occupationNotWorking_mobility_usa)

m.occupationNotWorking_food_germany <- glm(occupationGRP=="No paid working" ~ Food_, data= subset (df, country =="Germany" & group =="Food"), family=binomial())
summary(m.occupationNotWorking_food_germany)

m.occupationNotWorking_mobility_germany <- glm(occupationGRP=="No paid working" ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility" ), family=binomial())
summary(m.occupationNotWorking_mobility_germany)

m.occupationNotWorking_food_china <- glm(occupationGRP=="No paid working" ~ Food_, data= subset (df, country =="China" & group =="Food"), family=binomial())
summary(m.occupationNotWorking_food_china)

m.occupationNotWorking_mobility_china <- glm(occupationGRP=="No paid working" ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ), family=binomial())
summary(m.occupationNotWorking_mobility_china)

BToccupationNotWorking <- stargazer(m.occupationNotWorking_food_usa, m.occupationNotWorking_mobility_usa, m.occupationNotWorking_food_germany, m.occupationNotWorking_mobility_germany, m.occupationNotWorking_food_china, m.occupationNotWorking_mobility_china,
                                 type="html", 
                                 single.row = TRUE,
                                 title="Balance Test Not Working", 
                                 align=TRUE,
                                 star.char = c("+", "*", "**", "***"),
                                 star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                                 digits = 3,
                                 column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                                 object.names = FALSE,
                                 notes.append = FALSE, notes.align = "l",
                                 notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                                 out='BToccupationNotWorking.html')


m.pv_altr_group <- lm(pv_altr ~ group, df)
summary(m.pv_altr_group)

m.pv_altr_food_usa <- lm(pv_altr ~ Food_, data= subset (df, country =="USA" & group =="Food"))
summary(m.pv_altr_food_usa)

m.pv_altr_mobility_usa <- lm(pv_altr ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ))
summary(m.pv_altr_mobility_usa)

m.pv_altr_food_germany <- lm(pv_altr ~ Food_, data= subset (df, country =="Germany" & group =="Food"))
summary(m.pv_altr_food_germany)

m.pv_altr_mobility_germany <- lm(pv_altr ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility" ))
summary(m.pv_altr_mobility_germany)

m.pv_altr_food_china <- lm(pv_altr ~ Food_, data= subset (df, country =="China" & group =="Food"))
summary(m.pv_altr_food_china)

m.pv_altr_mobility_china <- lm(pv_altr ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ))
summary(m.pv_altr_mobility_china)

BTpv_altr <- stargazer(m.pv_altr_food_usa, m.pv_altr_mobility_usa, m.pv_altr_food_germany, m.pv_altr_mobility_germany, m.pv_altr_food_china, m.pv_altr_mobility_china,
                       type="html", 
                       single.row = TRUE,
                       title="Balance Test PV Altruistic", 
                       align=TRUE,
                       star.char = c("+", "*", "**", "***"),
                       star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                       digits = 3,
                       column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                       object.names = FALSE,
                       notes.append = FALSE, notes.align = "l",
                       notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                       out='BTpv_altr.html')

m.pv_ego_group <- lm(pv_ego ~ group, df)
summary(m.pv_ego_group)

m.pv_ego_food_usa <- lm(pv_ego ~ Food_, data= subset (df,country =="USA" & group =="Food"))
summary(m.pv_ego_food_usa)

m.pv_ego_mobility_usa <- lm(pv_ego ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ))
summary(m.pv_ego_mobility_usa)

m.pv_ego_food_germany <- lm(pv_ego ~ Food_, data= subset (df,country =="Germany" & group =="Food"))
summary(m.pv_ego_food_germany)

m.pv_ego_mobility_germany <- lm(pv_ego ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility" ))
summary(m.pv_ego_mobility_germany)

m.pv_ego_food_china <- lm(pv_ego ~ Food_, data= subset (df,country =="China" & group =="Food"))
summary(m.pv_ego_food_china)

m.pv_ego_mobility_china <- lm(pv_ego ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ))
summary(m.pv_ego_mobility_china)

BTpv_ego <- stargazer(m.pv_ego_food_usa, m.pv_ego_mobility_usa, m.pv_ego_food_germany, m.pv_ego_mobility_germany, m.pv_ego_food_china, m.pv_ego_mobility_china,
                      type="html", 
                      single.row = TRUE,
                      title="Balance Test PV Egoistic", 
                      align=TRUE,
                      star.char = c("+", "*", "**", "***"),
                      star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                      digits = 3,
                      column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                      object.names = FALSE,
                      notes.append = FALSE, notes.align = "l",
                      notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                      out='BTpv_ego.html')

m.pv_bios_group <- lm(pv_bios ~ group, df)
summary(m.pv_bios_group)

m.pv_bios_food_usa <- lm(pv_bios ~ Food_, data= subset (df, country =="USA" & group =="Food"))
summary(m.pv_bios_food_usa)

m.pv_bios_mobility_usa <- lm(pv_bios ~ Transport_, data= subset (df, country =="USA" & group =="Mobility" ))
summary(m.pv_bios_mobility_usa)

m.pv_bios_food_germany <- lm(pv_bios ~ Food_, data= subset (df, country =="Germany" & group =="Food"))
summary(m.pv_bios_food_germany)

m.pv_bios_mobility_germany <- lm(pv_bios ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility" ))
summary(m.pv_bios_mobility_germany)

m.pv_bios_food_china <- lm(pv_bios ~ Food_, data= subset (df, country =="China" & group =="Food"))
summary(m.pv_bios_food_china)

m.pv_bios_mobility_china <- lm(pv_bios ~ Transport_, data= subset (df, country =="China" & group =="Mobility" ))
summary(m.pv_bios_mobility_china)

BTpv_bios <- stargazer(m.pv_bios_food_usa, m.pv_bios_mobility_usa, m.pv_bios_food_germany, m.pv_bios_mobility_germany, m.pv_bios_food_china, m.pv_bios_mobility_china,
                       type="html", 
                       single.row = TRUE,
                       title="Balance Test PV Biospheric", 
                       align=TRUE,
                       star.char = c("+", "*", "**", "***"),
                       star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                       digits = 3,
                       column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                       object.names = FALSE,
                       notes.append = FALSE, notes.align = "l",
                       notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                       out='BTpv_bios.html')

m.pv_hedo_group <- lm(pv_hedo ~ group, df)
summary(m.pv_hedo_group)

m.pv_hedo_food_usa <- lm(pv_hedo ~ Food_, data= subset (df, country =="USA" & group =="Food"))
summary(m.pv_hedo_food_usa)

m.pv_hedo_mobility_usa <- lm(pv_hedo ~ Transport_, data= subset (df, country =="USA" & group =="Mobility"))
summary(m.pv_hedo_mobility_usa)

m.pv_hedo_food_germany <- lm(pv_hedo ~ Food_, data= subset (df, country =="Germany" & group =="Food"))
summary(m.pv_hedo_food_germany)

m.pv_hedo_mobility_germany <- lm(pv_hedo ~ Transport_, data= subset (df, country =="Germany" & group =="Mobility"))
summary(m.pv_hedo_mobility_germany)

m.pv_hedo_food_china <- lm(pv_hedo ~ Food_, data= subset (df, country =="China" & group =="Food"))
summary(m.pv_hedo_food_china)

m.pv_hedo_mobility_china <- lm(pv_hedo ~ Transport_, data= subset (df, country =="China" & group =="Mobility"))
summary(m.pv_hedo_mobility_china)

BTpv_hedo <- stargazer(m.pv_hedo_food_usa, m.pv_hedo_mobility_usa, m.pv_hedo_food_germany, m.pv_hedo_mobility_germany, m.pv_hedo_food_china, m.pv_hedo_mobility_china,
                       type="html", 
                       single.row = TRUE,
                       title="Balance Test PV Hedonistic", 
                       align=TRUE,
                       star.char = c("+", "*", "**", "***"),
                       star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                       digits = 3,
                       column.labels = c("USA Food", "USA Mobility", "Germany Food", "Germany Mobility", "China Food", "China Mobility"),
                       object.names = FALSE,
                       notes.append = FALSE, notes.align = "l",
                       notes = "Table entries are unstandardised regression coefficients from an OLS regression with estimated standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1",
                       out='BTpv_hedo.html')



#Prepare for LassoPlus Models---------

#Food dataset without NAs

id <- 1:nrow(df)
df_food_na <- data.frame(id)

df_food_na$policysupport_reduction_meatcon_z <- as.double(as.numeric(df$policysupport_reduction_meatcon_z))
df_food_na$slider_reduce_per_meatcon_z <- as.double(as.numeric(df$slider_reduce_per_meatcon_z))
df_food_na$slider_WTP_meat_z <- as.double(as.numeric(df$slider_WTP_meat_z))
df_food_na$concern_z <- as.double(as.numeric(df$concern_egotrophic_food_add_z))

df_food_na$pv_bios <- as.double(as.numeric(df$pv_bios))
df_food_na$pv_altr <- as.double(as.numeric(df$pv_altr))
df_food_na$pv_ego <- as.double(as.numeric(df$pv_ego))
df_food_na$pv_hedo <- as.double(as.numeric(df$pv_hedo))

df_food_na$gender <-as.double(df$gender)
df_food_na$age <-df$age
df_food_na$education <-ifelse(df$educationGRP=="Low", 1,
                                     ifelse(df$educationGRP=="Medium",2,3))
df_food_na$education <- as.double(df_food_na$education)
df_food_na$inc_n <-as.double(df$inc_n)
df_food_na$household_children <-as.double(df$household_children)
df_food_na$occupation <-ifelse(df$occupationGRP =="No paid working",0,1)
df_food_na$occupation <-as.double(df_food_na$occupation )

df_food_na$food_habits  <-as.double(as.numeric(as.factor(df$food_habits)))
df_food_na$meals_week <-as.double(as.numeric(df$meals_week))
df_food_na$amount_meat <-as.double(as.numeric(df$amount_meat))
df_food_na$dif_reduce <-as.double(as.numeric(df$dif_reduce))
df_food_na$quality_life <-as.double(as.numeric(df$quality_life))
df_food_na$cri_sust <-as.double(df$cri_sust)
df_food_na$cri_ego <-as.double(df$cri_ego)



df_food_na$Food_ <- as.factor(df$Food_)
df_food_na$treatment_food <- as.factor(df$treatment_food)

df_food_na$country <- as.factor(df$country)
df_food_na$country_n <- as.double(as.numeric(df_food_na$country))
df_food_na$group <- as.factor(df$group)

df_food_China_na <- subset(df_food_na,group=="Food" & country=="China")
df_food_China_na<- na.omit(df_food_China_na) 


#Political preference variables only asked in Germany and USA
df_food_na$government_int <-as.double(as.numeric(df$government_int))
df_food_na$left_right <-as.double(as.numeric(df$left_right))

#Urban/Rural variable only avaliable for Germany and USA
df_food_na$urbanRural  <-as.double(as.numeric(as.factor(df$urbanRural)))

df_food_USA_na <- subset(df_food_na,group=="Food" & country=="USA")
df_food_USA_na<- na.omit(df_food_USA_na) 
df_food_Germany_na <- subset(df_food_na,group=="Food" & country=="Germany")
df_food_Germany_na<- na.omit(df_food_Germany_na) 

#Transport dataset without NAs
id <- 1:nrow(df)
df_Mobility_na <- data.frame(id)

df_Mobility_na$policysupport_reduction_caruse_z <- as.double(as.numeric(df$policysupport_reduction_caruse_z))
df_Mobility_na$slider_reduce_per_caruse_z <- as.double(as.numeric(df$slider_reduce_per_caruse_z))
df_Mobility_na$slider_WTP_mobility_z <- as.double(as.numeric(df$slider_WTP_mobility_z))
df_Mobility_na$concern_z <- as.double(as.numeric(df$concern_egotrophic_Mobility_add_z))

df_Mobility_na$pv_bios <- as.double(as.numeric(df$pv_bios))
df_Mobility_na$pv_altr <- as.double(as.numeric(df$pv_altr))
df_Mobility_na$pv_ego <- as.double(as.numeric(df$pv_ego))
df_Mobility_na$pv_hedo <- as.double(as.numeric(df$pv_hedo))

df_Mobility_na$gender <-as.double(df$gender)
df_Mobility_na$age <-df$age
df_Mobility_na$education <-ifelse(df$educationGRP=="Low", 1,
                         ifelse(df$educationGRP=="Medium",2,3))
df_Mobility_na$education <- as.double(df_Mobility_na$education)
df_Mobility_na$inc_n <-as.double(df$inc_n)
df_Mobility_na$household_children <-as.double(df$household_children)
df_Mobility_na$occupation <-ifelse(df$occupationGRP =="No paid working",0,1)
df_Mobility_na$occupation <-as.double(df_Mobility_na$occupation )

df_Mobility_na$driving_habits_n   <-as.double(as.numeric(as.factor(df$driving_habits_n)))
df_Mobility_na$drivers_license  <-as.double(as.numeric(as.factor(df$drivers_license)))
df_Mobility_na$car_access  <-as.double(as.numeric(as.factor(df$car_access)))
df_Mobility_na$car_engine  <-as.double(as.numeric(as.factor(df$car_engine)))
df_Mobility_na$km_driven  <-as.double(as.numeric(df$km_driven ))
df_Mobility_na$crim_sust <-as.double(df$crim_sust)
df_Mobility_na$crim_ego <-as.double(df$crim_ego)

df_Mobility_na$dif_drivstop <-as.double(as.numeric(df$dif_drivstop))
df_Mobility_na$quality_life_m<-as.double(as.numeric(df$quality_life_m))

df_Mobility_na$Transport_ <- as.factor(df$Transport_)
df_Mobility_na$treatment_mobility <- as.factor(df$treatment_mobility)


df_Mobility_na$country <- as.factor(df$country)
df_Mobility_na$group <- as.factor(df$group)

df_Mobility_China_na <- subset(df_Mobility_na,group=="Mobility" & country=="China")
df_Mobility_China_na<- na.omit(df_Mobility_China_na) 

#Political preference variables only asked in Germany and USA
df_Mobility_na$government_int <-as.double(as.numeric(df$government_int))
df_Mobility_na$left_right <-as.double(as.numeric(df$left_right))

#Urban/Rural variable only avaliable for Germany and USA
df_Mobility_na$urbanRural  <-as.double(as.numeric(as.factor(df$urbanRural)))

df_Mobility_USA_na <- subset(df_Mobility_na,group=="Mobility" & country=="USA")
df_Mobility_USA_na<- na.omit(df_Mobility_USA_na) 
df_Mobility_Germany_na <- subset(df_Mobility_na,group=="Mobility" & country=="Germany")
df_Mobility_Germany_na<- na.omit(df_Mobility_Germany_na) 

#Estimate LassoPlus Models for Restricted Sample (Only individuals that eat meat/fish or drive a car themselves)--------------------------
#Please be aware that the estimation of the following LASSOplus regressions takes considerable amount of time

#Model 1 - DV - Concern about Impact of Meat/Fish Consumption/Car Use----------

#Food
out_new1_food_USA <- sparsereg::sparsereg(y = df_food_USA_na$concern_z, 
          X = as.matrix(data.frame(cri_sust = df_food_USA_na$cri_sust,cri_ego = df_food_USA_na$cri_ego,pv_ego= df_food_USA_na$pv_ego, pv_altr=df_food_USA_na$pv_altr,pv_bios= df_food_USA_na$pv_bios, pv_hedo=df_food_USA_na$pv_hedo,food_habits = df_food_USA_na$food_habits,meals_week = df_food_USA_na$meals_week,amount_meat = df_food_USA_na$amount_meat,dif_reduce = df_food_USA_na$dif_reduce,quality_life = df_food_USA_na$quality_life,occupation = df_food_USA_na$occupation,gender = df_food_USA_na$gender,age = df_food_USA_na$age,education = df_food_USA_na$education,inc_n = df_food_USA_na$inc_n,children = df_food_USA_na$household_children,urbanRural = df_food_USA_na$urbanRural, left_right = df_food_USA_na$left_right, government_int = df_food_USA_na$government_int)),
          treat = df_food_USA_na$Food_,
                            scale.type = "TTX",
          baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new1_food_USA)

sum1_food_USAa<- print(out_new1_food_USA)
sum1_food_USAb <- as.data.frame(sum1_food_USAa$statistics)
sum1_food_USAc <- as.data.frame(sum1_food_USAa$quantiles)
sum1_food_USA<-cbind (sum1_food_USAb,sum1_food_USAc)
sum1_food_USA$var <- rownames(sum1_food_USA)

sum1_food_USA<-sum1_food_USA %>%mutate_each(funs(round(., 3)), -var)
sum1_food_USA <- sum1_food_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum1_food_USA, file = "sum1_food_USA.csv", sep = ";", quote = FALSE, row.names = F)

sum1_food_USA <- sum1_food_USA [1:4,]
sum1_food_USA$country <- "USA"
sum1_food_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum1_food_USA$DV = factor("Concern about Meat/Fish Consumption")

out_new1_food_Germany <- sparsereg::sparsereg(y = df_food_Germany_na$concern_z, 
                                          X = as.matrix(data.frame(cri_sust = df_food_Germany_na$cri_sust,cri_ego = df_food_Germany_na$cri_ego,pv_ego= df_food_Germany_na$pv_ego, pv_altr=df_food_Germany_na$pv_altr,pv_bios= df_food_Germany_na$pv_bios, pv_hedo=df_food_Germany_na$pv_hedo,food_habits = df_food_Germany_na$food_habits,meals_week = df_food_Germany_na$meals_week,amount_meat = df_food_Germany_na$amount_meat,dif_reduce = df_food_Germany_na$dif_reduce,quality_life = df_food_Germany_na$quality_life,occupation = df_food_Germany_na$occupation,gender = df_food_Germany_na$gender,age = df_food_Germany_na$age,education = df_food_Germany_na$education,inc_n = df_food_Germany_na$inc_n,children = df_food_Germany_na$household_children,urbanRural = df_food_Germany_na$urbanRural, left_right = df_food_Germany_na$left_right, government_int = df_food_Germany_na$government_int)),
                                          treat = df_food_Germany_na$Food_,
                                          scale.type = "TTX",
                                          baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new1_food_Germany)

sum1_food_Germanya<- print(out_new1_food_Germany)
sum1_food_Germanyb <- as.data.frame(sum1_food_Germanya$statistics)
sum1_food_Germanyc <- as.data.frame(sum1_food_Germanya$quantiles)
sum1_food_Germany<-cbind (sum1_food_Germanyb,sum1_food_Germanyc)
sum1_food_Germany$var <- rownames(sum1_food_Germany)

sum1_food_Germany<-sum1_food_Germany %>%mutate_each(funs(round(., 3)), -var)
sum1_food_Germany <- sum1_food_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum1_food_Germany, file = "sum1_food_Germany.csv", sep = ";", quote = FALSE, row.names = F)

sum1_food_Germany <- sum1_food_Germany [1:4,]
sum1_food_Germany$country <- "Germany"
sum1_food_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum1_food_Germany$DV = factor("Concern about Meat/Fish Consumption")

out_new1_food_China <- sparsereg::sparsereg(y = df_food_China_na$concern_z, 
                                              X = as.matrix(data.frame(cri_sust = df_food_China_na$cri_sust,cri_ego = df_food_China_na$cri_ego,pv_ego= df_food_China_na$pv_ego, pv_altr=df_food_China_na$pv_altr,pv_bios= df_food_China_na$pv_bios, pv_hedo=df_food_China_na$pv_hedo,food_habits = df_food_China_na$food_habits,meals_week = df_food_China_na$meals_week,amount_meat = df_food_China_na$amount_meat,dif_reduce = df_food_China_na$dif_reduce,quality_life = df_food_China_na$quality_life,occupation = df_food_China_na$occupation,gender = df_food_China_na$gender,age = df_food_China_na$age,education = df_food_China_na$education,inc_n = df_food_China_na$inc_n,children = df_food_China_na$household_children)),
                                              treat = df_food_China_na$Food_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new1_food_China)

sum1_food_Chinaa<- print(out_new1_food_China)
sum1_food_Chinab <- as.data.frame(sum1_food_Chinaa$statistics)
sum1_food_Chinac <- as.data.frame(sum1_food_Chinaa$quantiles)
sum1_food_China<-cbind (sum1_food_Chinab,sum1_food_Chinac)
sum1_food_China$var <- rownames(sum1_food_China)

sum1_food_China<-sum1_food_China %>%mutate_each(funs(round(., 3)), -var)
sum1_food_China <- sum1_food_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum1_food_China, file = "sum1_food_China.csv", sep = ";", quote = FALSE, row.names = F)

sum1_food_China <- sum1_food_China [1:4,]
sum1_food_China$country <- "China"
sum1_food_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum1_food_China$DV = factor("Concern about Meat/Fish Consumption")

#Mobility

out_new1_Mobility_USA <- sparsereg::sparsereg(y = df_Mobility_USA_na$concern_z, 
                                             X = as.matrix(data.frame(crim_sust = df_Mobility_USA_na$crim_sust,crim_ego = df_Mobility_USA_na$crim_ego,pv_ego= df_Mobility_USA_na$pv_ego, pv_altr=df_Mobility_USA_na$pv_altr,pv_bios= df_Mobility_USA_na$pv_bios, pv_hedo=df_Mobility_USA_na$pv_hedo,driving_habits = df_Mobility_USA_na$driving_habits,drivers_license=df_Mobility_USA_na$drivers_license,car_access= df_Mobility_USA_na$car_access,car_engine=df_Mobility_USA_na$car_engine,  km_driven  = df_Mobility_USA_na$km_driven,dif_drivstop=df_Mobility_USA_na$dif_drivstop, quality_life_m=df_Mobility_USA_na$quality_life_m, occupation = df_Mobility_USA_na$occupation,gender = df_Mobility_USA_na$gender,age = df_Mobility_USA_na$age,education = df_Mobility_USA_na$education,inc_n = df_Mobility_USA_na$inc_n,children = df_Mobility_USA_na$household_children,urbanRural = df_Mobility_USA_na$urbanRural, left_right = df_Mobility_USA_na$left_right, government_int = df_Mobility_USA_na$government_int)),
                                             treat = df_Mobility_USA_na$Transport_,
                                             scale.type = "TTX",
                                             baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new1_Mobility_USA)

sum1_Mobility_USAa<- print(out_new1_Mobility_USA)
sum1_Mobility_USAb <- as.data.frame(sum1_Mobility_USAa$statistics)
sum1_Mobility_USAc <- as.data.frame(sum1_Mobility_USAa$quantiles)
sum1_Mobility_USA<-cbind (sum1_Mobility_USAb,sum1_Mobility_USAc)
sum1_Mobility_USA$var <- rownames(sum1_Mobility_USA)

sum1_Mobility_USA<-sum1_Mobility_USA %>%mutate_each(funs(round(., 3)), -var)
sum1_Mobility_USA <- sum1_Mobility_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum1_Mobility_USA, file = "sum1_Mobility_USA.csv", sep = ";", quote = FALSE, row.names = F)


sum1_Mobility_USA <- sum1_Mobility_USA [1:4,]
sum1_Mobility_USA$country <- "USA"
sum1_Mobility_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum1_Mobility_USA$DV = factor("Concern about Car Use")


out_new1_Mobility_Germany <- sparsereg::sparsereg(y = df_Mobility_Germany_na$concern_z, 
                                              X = as.matrix(data.frame(crim_sust = df_Mobility_Germany_na$crim_sust,crim_ego = df_Mobility_Germany_na$crim_ego,pv_ego= df_Mobility_Germany_na$pv_ego, pv_altr=df_Mobility_Germany_na$pv_altr,pv_bios= df_Mobility_Germany_na$pv_bios, pv_hedo=df_Mobility_Germany_na$pv_hedo,driving_habits = df_Mobility_Germany_na$driving_habits,drivers_license=df_Mobility_Germany_na$drivers_license,car_access= df_Mobility_Germany_na$car_access,car_engine=df_Mobility_Germany_na$car_engine,  km_driven  = df_Mobility_Germany_na$km_driven,dif_drivstop=df_Mobility_Germany_na$dif_drivstop, quality_life_m=df_Mobility_Germany_na$quality_life_m, occupation = df_Mobility_Germany_na$occupation,gender = df_Mobility_Germany_na$gender,age = df_Mobility_Germany_na$age,education = df_Mobility_Germany_na$education,inc_n = df_Mobility_Germany_na$inc_n,children = df_Mobility_Germany_na$household_children,urbanRural = df_Mobility_Germany_na$urbanRural, left_right = df_Mobility_Germany_na$left_right, government_int = df_Mobility_Germany_na$government_int)),
                                              treat = df_Mobility_Germany_na$Transport_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new1_Mobility_Germany)

sum1_Mobility_Germanya<- print(out_new1_Mobility_Germany)
sum1_Mobility_Germanyb <- as.data.frame(sum1_Mobility_Germanya$statistics)
sum1_Mobility_Germanyc <- as.data.frame(sum1_Mobility_Germanya$quantiles)
sum1_Mobility_Germany<-cbind (sum1_Mobility_Germanyb,sum1_Mobility_Germanyc)
sum1_Mobility_Germany$var <- rownames(sum1_Mobility_Germany)

sum1_Mobility_Germany<-sum1_Mobility_Germany %>%mutate_each(funs(round(., 3)), -var)
sum1_Mobility_Germany <- sum1_Mobility_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum1_Mobility_Germany, file = "sum1_Mobility_Germany.csv", sep = ";", quote = FALSE, row.names = F)

sum1_Mobility_Germany <- sum1_Mobility_Germany [1:4,]
sum1_Mobility_Germany$country <- "Germany"
sum1_Mobility_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum1_Mobility_Germany$DV = factor("Concern about Car Use")

out_new1_Mobility_China <- sparsereg::sparsereg(y = df_Mobility_China_na$concern_z, 
                                              X = as.matrix(data.frame(crim_sust = df_Mobility_China_na$crim_sust,crim_ego = df_Mobility_China_na$crim_ego,pv_ego= df_Mobility_China_na$pv_ego, pv_altr=df_Mobility_China_na$pv_altr,pv_bios= df_Mobility_China_na$pv_bios, pv_hedo=df_Mobility_China_na$pv_hedo,driving_habits = df_Mobility_China_na$driving_habits,drivers_license=df_Mobility_China_na$drivers_license,car_access= df_Mobility_China_na$car_access,car_engine=df_Mobility_China_na$car_engine,  km_driven  = df_Mobility_China_na$km_driven,dif_drivstop=df_Mobility_China_na$dif_drivstop, quality_life_m=df_Mobility_China_na$quality_life_m, occupation = df_Mobility_China_na$occupation,gender = df_Mobility_China_na$gender,age = df_Mobility_China_na$age,education = df_Mobility_China_na$education,inc_n = df_Mobility_China_na$inc_n,children = df_Mobility_China_na$household_children)),
                                              treat = df_Mobility_China_na$Transport_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new1_Mobility_China)

sum1_Mobility_Chinaa<- print(out_new1_Mobility_China)
sum1_Mobility_Chinab <- as.data.frame(sum1_Mobility_Chinaa$statistics)
sum1_Mobility_Chinac <- as.data.frame(sum1_Mobility_Chinaa$quantiles)
sum1_Mobility_China<-cbind (sum1_Mobility_Chinab,sum1_Mobility_Chinac)
sum1_Mobility_China$var <- rownames(sum1_Mobility_China)

sum1_Mobility_China<-sum1_Mobility_China %>%mutate_each(funs(round(., 3)), -var)
sum1_Mobility_China <- sum1_Mobility_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum1_Mobility_China, file = "sum1_Mobility_China.csv", sep = ";", quote = FALSE, row.names = F)

sum1_Mobility_China <- sum1_Mobility_China [1:4,]
sum1_Mobility_China$country <- "China"
sum1_Mobility_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum1_Mobility_China$DV = factor("Concern about Car Use")


sum1_df <- rbind(sum1_Mobility_China,sum1_Mobility_Germany,sum1_Mobility_USA,sum1_food_China,sum1_food_Germany, sum1_food_USA)


sum1_dftable <- cbind(sum1_Mobility_China,sum1_Mobility_Germany,sum1_Mobility_USA,sum1_food_China,sum1_food_Germany, sum1_food_USA)

#Model 2: DV - Q14.3.	SA PER ITEM/ Policy Support Reduction----------

summary(df_food_USA$policysupport_reduction_meatcon [df_food_USA$treatment_food=="Control"])



#Food
out_new2_food_USA <- sparsereg::sparsereg(y = df_food_USA_na$policysupport_reduction_meatcon_z, 
                                          X = as.matrix(data.frame(cri_sust = df_food_USA_na$cri_sust,cri_ego = df_food_USA_na$cri_ego,pv_ego= df_food_USA_na$pv_ego, pv_altr=df_food_USA_na$pv_altr,pv_bios= df_food_USA_na$pv_bios, pv_hedo=df_food_USA_na$pv_hedo,food_habits = df_food_USA_na$food_habits,meals_week = df_food_USA_na$meals_week,amount_meat = df_food_USA_na$amount_meat,dif_reduce = df_food_USA_na$dif_reduce,quality_life = df_food_USA_na$quality_life,occupation = df_food_USA_na$occupation,gender = df_food_USA_na$gender,age = df_food_USA_na$age,education = df_food_USA_na$education,inc_n = df_food_USA_na$inc_n,children = df_food_USA_na$household_children,urbanRural = df_food_USA_na$urbanRural, left_right = df_food_USA_na$left_right, government_int = df_food_USA_na$government_int)),
                                          treat = df_food_USA_na$Food_,
                                          scale.type = "TTX",
                                          baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new2_food_USA)

sum2_food_USAa<- print(out_new2_food_USA)
sum2_food_USAb <- as.data.frame(sum2_food_USAa$statistics)
sum2_food_USAc <- as.data.frame(sum2_food_USAa$quantiles)
sum2_food_USA<-cbind (sum2_food_USAb,sum2_food_USAc)
sum2_food_USA$var <- rownames(sum2_food_USA)

sum2_food_USA<-sum2_food_USA %>%mutate_each(funs(round(., 3)), -var)
sum2_food_USA <- sum2_food_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum2_food_USA, file = "sum2_food_USA.csv", sep = ";", quote = FALSE, row.names = F)


sum2_food_USA <- sum2_food_USA [1:4,]
sum2_food_USA$country <- "USA"
sum2_food_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum2_food_USA$DV = factor("Policy Support to Reduce Meat/Fish Consumption")


out_new2_food_Germany <- sparsereg::sparsereg(y = df_food_Germany_na$policysupport_reduction_meatcon_z, 
                                              X = as.matrix(data.frame(cri_sust = df_food_Germany_na$cri_sust,cri_ego = df_food_Germany_na$cri_ego,pv_ego= df_food_Germany_na$pv_ego, pv_altr=df_food_Germany_na$pv_altr,pv_bios= df_food_Germany_na$pv_bios, pv_hedo=df_food_Germany_na$pv_hedo,food_habits = df_food_Germany_na$food_habits,meals_week = df_food_Germany_na$meals_week,amount_meat = df_food_Germany_na$amount_meat,dif_reduce = df_food_Germany_na$dif_reduce,quality_life = df_food_Germany_na$quality_life,occupation = df_food_Germany_na$occupation,gender = df_food_Germany_na$gender,age = df_food_Germany_na$age,education = df_food_Germany_na$education,inc_n = df_food_Germany_na$inc_n,children = df_food_Germany_na$household_children,urbanRural = df_food_Germany_na$urbanRural, left_right = df_food_Germany_na$left_right, government_int = df_food_Germany_na$government_int)),
                                              treat = df_food_Germany_na$Food_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new2_food_Germany)


sum2_food_Germanya<- print(out_new2_food_Germany)
sum2_food_Germanyb <- as.data.frame(sum2_food_Germanya$statistics)
sum2_food_Germanyc <- as.data.frame(sum2_food_Germanya$quantiles)
sum2_food_Germany<-cbind (sum2_food_Germanyb,sum2_food_Germanyc)
sum2_food_Germany$var <- rownames(sum2_food_Germany)

sum2_food_Germany<-sum2_food_Germany %>%mutate_each(funs(round(., 3)), -var)
sum2_food_Germany <- sum2_food_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum2_food_Germany, file = "sum2_food_Germany.csv", sep = ";", quote = FALSE, row.names = F)

sum2_food_Germany <- sum2_food_Germany [1:4,]
sum2_food_Germany$country <- "Germany"
sum2_food_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum2_food_Germany$DV = factor("Policy Support to Reduce Meat/Fish Consumption")



out_new2_food_China <- sparsereg::sparsereg(y = df_food_China_na$policysupport_reduction_meatcon_z, 
                                            X = as.matrix(data.frame(cri_sust = df_food_China_na$cri_sust,cri_ego = df_food_China_na$cri_ego,pv_ego= df_food_China_na$pv_ego, pv_altr=df_food_China_na$pv_altr,pv_bios= df_food_China_na$pv_bios, pv_hedo=df_food_China_na$pv_hedo,food_habits = df_food_China_na$food_habits,meals_week = df_food_China_na$meals_week,amount_meat = df_food_China_na$amount_meat,dif_reduce = df_food_China_na$dif_reduce,quality_life = df_food_China_na$quality_life,occupation = df_food_China_na$occupation,gender = df_food_China_na$gender,age = df_food_China_na$age,education = df_food_China_na$education,inc_n = df_food_China_na$inc_n,children = df_food_China_na$household_children)),
                                            treat = df_food_China_na$Food_,
                                            scale.type = "TTX",
                                            baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new2_food_China)


sum2_food_Chinaa<- print(out_new2_food_China)
sum2_food_Chinab <- as.data.frame(sum2_food_Chinaa$statistics)
sum2_food_Chinac <- as.data.frame(sum2_food_Chinaa$quantiles)
sum2_food_China<-cbind (sum2_food_Chinab,sum2_food_Chinac)
sum2_food_China$var <- rownames(sum2_food_China)

sum2_food_China<-sum2_food_China %>%mutate_each(funs(round(., 3)), -var)
sum2_food_China <- sum2_food_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum2_food_China, file = "sum2_food_China.csv", sep = ";", quote = FALSE, row.names = F)

sum2_food_China <- sum2_food_China [1:4,]
sum2_food_China$country <- "China"
sum2_food_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum2_food_China$DV = factor("Policy Support to Reduce Meat/Fish Consumption")

#Mobility
out_new2_Mobility_USA <- sparsereg::sparsereg(y = df_Mobility_USA_na$policysupport_reduction_caruse_z, 
                                              X = as.matrix(data.frame(crim_sust = df_Mobility_USA_na$crim_sust,crim_ego = df_Mobility_USA_na$crim_ego,pv_ego= df_Mobility_USA_na$pv_ego, pv_altr=df_Mobility_USA_na$pv_altr,pv_bios= df_Mobility_USA_na$pv_bios, pv_hedo=df_Mobility_USA_na$pv_hedo,driving_habits = df_Mobility_USA_na$driving_habits,drivers_license=df_Mobility_USA_na$drivers_license,car_access= df_Mobility_USA_na$car_access,car_engine=df_Mobility_USA_na$car_engine,  km_driven  = df_Mobility_USA_na$km_driven,dif_drivstop=df_Mobility_USA_na$dif_drivstop, quality_life_m=df_Mobility_USA_na$quality_life_m, occupation = df_Mobility_USA_na$occupation,gender = df_Mobility_USA_na$gender,age = df_Mobility_USA_na$age,education = df_Mobility_USA_na$education,inc_n = df_Mobility_USA_na$inc_n,children = df_Mobility_USA_na$household_children,urbanRural = df_Mobility_USA_na$urbanRural, left_right = df_Mobility_USA_na$left_right, government_int = df_Mobility_USA_na$government_int)),
                                              treat = df_Mobility_USA_na$Transport_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new2_Mobility_USA)

sum2_Mobility_USAa<- print(out_new2_Mobility_USA)
sum2_Mobility_USAb <- as.data.frame(sum2_Mobility_USAa$statistics)
sum2_Mobility_USAc <- as.data.frame(sum2_Mobility_USAa$quantiles)
sum2_Mobility_USA<-cbind (sum2_Mobility_USAb,sum2_Mobility_USAc)
sum2_Mobility_USA$var <- rownames(sum2_Mobility_USA)

sum2_Mobility_USA<-sum2_Mobility_USA %>%mutate_each(funs(round(., 3)), -var)
sum2_Mobility_USA <- sum2_Mobility_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum2_Mobility_USA, file = "sum2_Mobility_USA.csv", sep = ";", quote = FALSE, row.names = F)

sum2_Mobility_USA <- sum2_Mobility_USA [1:4,]
sum2_Mobility_USA$country <- "USA"
sum2_Mobility_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum2_Mobility_USA$DV = factor("Policy Support to Reduce Car Use")


out_new2_Mobility_Germany <- sparsereg::sparsereg(y = df_Mobility_Germany_na$policysupport_reduction_caruse_z, 
                                                  X = as.matrix(data.frame(crim_sust = df_Mobility_Germany_na$crim_sust,crim_ego = df_Mobility_Germany_na$crim_ego,pv_ego= df_Mobility_Germany_na$pv_ego, pv_altr=df_Mobility_Germany_na$pv_altr,pv_bios= df_Mobility_Germany_na$pv_bios, pv_hedo=df_Mobility_Germany_na$pv_hedo,driving_habits = df_Mobility_Germany_na$driving_habits,drivers_license=df_Mobility_Germany_na$drivers_license,car_access= df_Mobility_Germany_na$car_access,car_engine=df_Mobility_Germany_na$car_engine,  km_driven  = df_Mobility_Germany_na$km_driven,dif_drivstop=df_Mobility_Germany_na$dif_drivstop, quality_life_m=df_Mobility_Germany_na$quality_life_m, occupation = df_Mobility_Germany_na$occupation,gender = df_Mobility_Germany_na$gender,age = df_Mobility_Germany_na$age,education = df_Mobility_Germany_na$education,inc_n = df_Mobility_Germany_na$inc_n,children = df_Mobility_Germany_na$household_children,urbanRural = df_Mobility_Germany_na$urbanRural, left_right = df_Mobility_Germany_na$left_right, government_int = df_Mobility_Germany_na$government_int)),
                                                  treat = df_Mobility_Germany_na$Transport_,
                                                  scale.type = "TTX",
                                                  baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new2_Mobility_Germany)

sum2_Mobility_Germanya<- print(out_new2_Mobility_Germany)
sum2_Mobility_Germanyb <- as.data.frame(sum2_Mobility_Germanya$statistics)
sum2_Mobility_Germanyc <- as.data.frame(sum2_Mobility_Germanya$quantiles)
sum2_Mobility_Germany<-cbind (sum2_Mobility_Germanyb,sum2_Mobility_Germanyc)
sum2_Mobility_Germany$var <- rownames(sum2_Mobility_Germany)

sum2_Mobility_Germany<-sum2_Mobility_Germany %>%mutate_each(funs(round(., 3)), -var)
sum2_Mobility_Germany <- sum2_Mobility_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum2_Mobility_Germany, file = "sum2_Mobility_Germany.csv", sep = ";", quote = FALSE, row.names = F)

sum2_Mobility_Germany <- sum2_Mobility_Germany [1:4,]
sum2_Mobility_Germany$country <- "Germany"
sum2_Mobility_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum2_Mobility_Germany$DV = factor("Policy Support to Reduce Car Use")

out_new2_Mobility_China <- sparsereg::sparsereg(y = df_Mobility_China_na$policysupport_reduction_caruse_z, 
                                                X = as.matrix(data.frame(crim_sust = df_Mobility_China_na$crim_sust,crim_ego = df_Mobility_China_na$crim_ego,pv_ego= df_Mobility_China_na$pv_ego, pv_altr=df_Mobility_China_na$pv_altr,pv_bios= df_Mobility_China_na$pv_bios, pv_hedo=df_Mobility_China_na$pv_hedo,driving_habits = df_Mobility_China_na$driving_habits,drivers_license=df_Mobility_China_na$drivers_license,car_access= df_Mobility_China_na$car_access,car_engine=df_Mobility_China_na$car_engine,  km_driven  = df_Mobility_China_na$km_driven,dif_drivstop=df_Mobility_China_na$dif_drivstop, quality_life_m=df_Mobility_China_na$quality_life_m, occupation = df_Mobility_China_na$occupation,gender = df_Mobility_China_na$gender,age = df_Mobility_China_na$age,education = df_Mobility_China_na$education,inc_n = df_Mobility_China_na$inc_n,children = df_Mobility_China_na$household_children)),
                                                treat = df_Mobility_China_na$Transport_,
                                                scale.type = "TTX",
                                                baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new2_Mobility_China)

sum2_Mobility_Chinaa<- print(out_new2_Mobility_China)
sum2_Mobility_Chinab <- as.data.frame(sum2_Mobility_Chinaa$statistics)
sum2_Mobility_Chinac <- as.data.frame(sum2_Mobility_Chinaa$quantiles)
sum2_Mobility_China<-cbind (sum2_Mobility_Chinab,sum2_Mobility_Chinac)
sum2_Mobility_China$var <- rownames(sum2_Mobility_China)

sum2_Mobility_China<-sum2_Mobility_China %>%mutate_each(funs(round(., 3)), -var)
sum2_Mobility_China <- sum2_Mobility_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum2_Mobility_China, file = "sum2_Mobility_China.csv", sep = ";", quote = FALSE, row.names = F)


sum2_Mobility_China <- sum2_Mobility_China [1:4,]
sum2_Mobility_China$country <- "China"
sum2_Mobility_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum2_Mobility_China$DV = factor("Policy Support to Reduce Car Use")


sum2_df <- rbind(sum2_Mobility_China,sum2_Mobility_Germany,sum2_Mobility_USA,sum2_food_China,sum2_food_Germany, sum2_food_USA)


#Model 3: DV - Q14.1 Slider/Personal WTP--------
#Food
out_new3_food_USA <- sparsereg::sparsereg(y = df_food_USA_na$slider_WTP_meat_z, 
                                          X = as.matrix(data.frame(cri_sust = df_food_USA_na$cri_sust,cri_ego = df_food_USA_na$cri_ego,pv_ego= df_food_USA_na$pv_ego, pv_altr=df_food_USA_na$pv_altr,pv_bios= df_food_USA_na$pv_bios, pv_hedo=df_food_USA_na$pv_hedo,food_habits = df_food_USA_na$food_habits,meals_week = df_food_USA_na$meals_week,amount_meat = df_food_USA_na$amount_meat,dif_reduce = df_food_USA_na$dif_reduce,quality_life = df_food_USA_na$quality_life,occupation = df_food_USA_na$occupation,gender = df_food_USA_na$gender,age = df_food_USA_na$age,education = df_food_USA_na$education,inc_n = df_food_USA_na$inc_n,children = df_food_USA_na$household_children,urbanRural = df_food_USA_na$urbanRural, left_right = df_food_USA_na$left_right, government_int = df_food_USA_na$government_int)),
                                          treat = df_food_USA_na$Food_,
                                          scale.type = "TTX",
                                          baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new3_food_USA)

sum3_food_USAa<- print(out_new3_food_USA)
sum3_food_USAb <- as.data.frame(sum3_food_USAa$statistics)
sum3_food_USAc <- as.data.frame(sum3_food_USAa$quantiles)
sum3_food_USA<-cbind (sum3_food_USAb,sum3_food_USAc)
sum3_food_USA$var <- rownames(sum3_food_USA)

sum3_food_USA<-sum3_food_USA %>%mutate_each(funs(round(., 3)), -var)
sum3_food_USA <- sum3_food_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum3_food_USA, file = "sum3_food_USA.csv", sep = ";", quote = FALSE, row.names = F)


sum3_food_USA <- sum3_food_USA [1:4,]
sum3_food_USA$country <- "USA"
sum3_food_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum3_food_USA$DV = factor("Willingness to Pay more for Meat/Fish")


out_new3_food_Germany <- sparsereg::sparsereg(y = df_food_Germany_na$slider_WTP_meat_z, 
                                              X = as.matrix(data.frame(cri_sust = df_food_Germany_na$cri_sust,cri_ego = df_food_Germany_na$cri_ego,pv_ego= df_food_Germany_na$pv_ego, pv_altr=df_food_Germany_na$pv_altr,pv_bios= df_food_Germany_na$pv_bios, pv_hedo=df_food_Germany_na$pv_hedo,food_habits = df_food_Germany_na$food_habits,meals_week = df_food_Germany_na$meals_week,amount_meat = df_food_Germany_na$amount_meat,dif_reduce = df_food_Germany_na$dif_reduce,quality_life = df_food_Germany_na$quality_life,occupation = df_food_Germany_na$occupation,gender = df_food_Germany_na$gender,age = df_food_Germany_na$age,education = df_food_Germany_na$education,inc_n = df_food_Germany_na$inc_n,children = df_food_Germany_na$household_children,urbanRural = df_food_Germany_na$urbanRural, left_right = df_food_Germany_na$left_right, government_int = df_food_Germany_na$government_int)),
                                              treat = df_food_Germany_na$Food_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new3_food_Germany)


sum3_food_Germanya<- print(out_new3_food_Germany)
sum3_food_Germanyb <- as.data.frame(sum3_food_Germanya$statistics)
sum3_food_Germanyc <- as.data.frame(sum3_food_Germanya$quantiles)
sum3_food_Germany<-cbind (sum3_food_Germanyb,sum3_food_Germanyc)
sum3_food_Germany$var <- rownames(sum3_food_Germany)

sum3_food_Germany<-sum3_food_Germany %>%mutate_each(funs(round(., 3)), -var)
sum3_food_Germany <- sum3_food_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum3_food_Germany, file = "sum3_food_Germany.csv", sep = ";", quote = FALSE, row.names = F)


sum3_food_Germany <- sum3_food_Germany [1:4,]
sum3_food_Germany$country <- "Germany"
sum3_food_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum3_food_Germany$DV = factor("Willingness to Pay more for Meat/Fish")



out_new3_food_China <- sparsereg::sparsereg(y = df_food_China_na$slider_WTP_meat_z, 
                                            X = as.matrix(data.frame(cri_sust = df_food_China_na$cri_sust,cri_ego = df_food_China_na$cri_ego,pv_ego= df_food_China_na$pv_ego, pv_altr=df_food_China_na$pv_altr,pv_bios= df_food_China_na$pv_bios, pv_hedo=df_food_China_na$pv_hedo,food_habits = df_food_China_na$food_habits,meals_week = df_food_China_na$meals_week,amount_meat = df_food_China_na$amount_meat,dif_reduce = df_food_China_na$dif_reduce,quality_life = df_food_China_na$quality_life,occupation = df_food_China_na$occupation,gender = df_food_China_na$gender,age = df_food_China_na$age,education = df_food_China_na$education,inc_n = df_food_China_na$inc_n,children = df_food_China_na$household_children)),
                                            treat = df_food_China_na$Food_,
                                            scale.type = "TTX",
                                            baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new3_food_China)


sum3_food_Chinaa<- print(out_new3_food_China)
sum3_food_Chinab <- as.data.frame(sum3_food_Chinaa$statistics)
sum3_food_Chinac <- as.data.frame(sum3_food_Chinaa$quantiles)
sum3_food_China<-cbind (sum3_food_Chinab,sum3_food_Chinac)
sum3_food_China$var <- rownames(sum3_food_China)

sum3_food_China<-sum3_food_China %>%mutate_each(funs(round(., 3)), -var)
sum3_food_China <- sum3_food_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum3_food_China, file = "sum3_food_China.csv", sep = ";", quote = FALSE, row.names = F)


sum3_food_China <- sum3_food_China [1:4,]
sum3_food_China$country <- "China"
sum3_food_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum3_food_China$DV = factor("Willingness to Pay more for Meat/Fish")

#Mobility
out_new3_Mobility_USA <- sparsereg::sparsereg(y = df_Mobility_USA_na$slider_WTP_mobility_z, 
                                              X = as.matrix(data.frame(crim_sust = df_Mobility_USA_na$crim_sust,crim_ego = df_Mobility_USA_na$crim_ego,pv_ego= df_Mobility_USA_na$pv_ego, pv_altr=df_Mobility_USA_na$pv_altr,pv_bios= df_Mobility_USA_na$pv_bios, pv_hedo=df_Mobility_USA_na$pv_hedo,driving_habits = df_Mobility_USA_na$driving_habits,drivers_license=df_Mobility_USA_na$drivers_license,car_access= df_Mobility_USA_na$car_access,car_engine=df_Mobility_USA_na$car_engine,  km_driven  = df_Mobility_USA_na$km_driven,dif_drivstop=df_Mobility_USA_na$dif_drivstop, quality_life_m=df_Mobility_USA_na$quality_life_m, occupation = df_Mobility_USA_na$occupation,gender = df_Mobility_USA_na$gender,age = df_Mobility_USA_na$age,education = df_Mobility_USA_na$education,inc_n = df_Mobility_USA_na$inc_n,children = df_Mobility_USA_na$household_children,urbanRural = df_Mobility_USA_na$urbanRural, left_right = df_Mobility_USA_na$left_right, government_int = df_Mobility_USA_na$government_int)),
                                              treat = df_Mobility_USA_na$Transport_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new3_Mobility_USA)

sum3_Mobility_USAa<- print(out_new3_Mobility_USA)
sum3_Mobility_USAb <- as.data.frame(sum3_Mobility_USAa$statistics)
sum3_Mobility_USAc <- as.data.frame(sum3_Mobility_USAa$quantiles)
sum3_Mobility_USA<-cbind (sum3_Mobility_USAb,sum3_Mobility_USAc)
sum3_Mobility_USA$var <- rownames(sum3_Mobility_USA)

sum3_Mobility_USA<-sum3_Mobility_USA %>%mutate_each(funs(round(., 3)), -var)
sum3_Mobility_USA <- sum3_Mobility_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum3_Mobility_USA, file = "sum3_Mobility_USA.csv", sep = ";", quote = FALSE, row.names = F)

sum3_Mobility_USA <- sum3_Mobility_USA [1:4,]
sum3_Mobility_USA$country <- "USA"
sum3_Mobility_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum3_Mobility_USA$DV = factor("Willingness to Pay more for Motor Fuels")


out_new3_Mobility_Germany <- sparsereg::sparsereg(y = df_Mobility_Germany_na$slider_WTP_mobility_z, 
                                                  X = as.matrix(data.frame(crim_sust = df_Mobility_Germany_na$crim_sust,crim_ego = df_Mobility_Germany_na$crim_ego,pv_ego= df_Mobility_Germany_na$pv_ego, pv_altr=df_Mobility_Germany_na$pv_altr,pv_bios= df_Mobility_Germany_na$pv_bios, pv_hedo=df_Mobility_Germany_na$pv_hedo,driving_habits = df_Mobility_Germany_na$driving_habits,drivers_license=df_Mobility_Germany_na$drivers_license,car_access= df_Mobility_Germany_na$car_access,car_engine=df_Mobility_Germany_na$car_engine,  km_driven  = df_Mobility_Germany_na$km_driven,dif_drivstop=df_Mobility_Germany_na$dif_drivstop, quality_life_m=df_Mobility_Germany_na$quality_life_m, occupation = df_Mobility_Germany_na$occupation,gender = df_Mobility_Germany_na$gender,age = df_Mobility_Germany_na$age,education = df_Mobility_Germany_na$education,inc_n = df_Mobility_Germany_na$inc_n,children = df_Mobility_Germany_na$household_children,urbanRural = df_Mobility_Germany_na$urbanRural, left_right = df_Mobility_Germany_na$left_right, government_int = df_Mobility_Germany_na$government_int)),
                                                  treat = df_Mobility_Germany_na$Transport_,
                                                  scale.type = "TTX",
                                                  baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new3_Mobility_Germany)

sum3_Mobility_Germanya<- print(out_new3_Mobility_Germany)
sum3_Mobility_Germanyb <- as.data.frame(sum3_Mobility_Germanya$statistics)
sum3_Mobility_Germanyc <- as.data.frame(sum3_Mobility_Germanya$quantiles)
sum3_Mobility_Germany<-cbind (sum3_Mobility_Germanyb,sum3_Mobility_Germanyc)
sum3_Mobility_Germany$var <- rownames(sum3_Mobility_Germany)

sum3_Mobility_Germany<-sum3_Mobility_Germany %>%mutate_each(funs(round(., 3)), -var)
sum3_Mobility_Germany <- sum3_Mobility_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum3_Mobility_Germany, file = "sum3_Mobility_Germany.csv", sep = ";", quote = FALSE, row.names = F)

sum3_Mobility_Germany <- sum3_Mobility_Germany [1:4,]
sum3_Mobility_Germany$country <- "Germany"
sum3_Mobility_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum3_Mobility_Germany$DV = factor("Willingness to Pay more for Motor Fuels")

out_new3_Mobility_China <- sparsereg::sparsereg(y = df_Mobility_China_na$slider_WTP_mobility_z, 
                                                X = as.matrix(data.frame(crim_sust = df_Mobility_China_na$crim_sust,crim_ego = df_Mobility_China_na$crim_ego,pv_ego= df_Mobility_China_na$pv_ego, pv_altr=df_Mobility_China_na$pv_altr,pv_bios= df_Mobility_China_na$pv_bios, pv_hedo=df_Mobility_China_na$pv_hedo,driving_habits = df_Mobility_China_na$driving_habits,drivers_license=df_Mobility_China_na$drivers_license,car_access= df_Mobility_China_na$car_access,car_engine=df_Mobility_China_na$car_engine,  km_driven  = df_Mobility_China_na$km_driven,dif_drivstop=df_Mobility_China_na$dif_drivstop, quality_life_m=df_Mobility_China_na$quality_life_m, occupation = df_Mobility_China_na$occupation,gender = df_Mobility_China_na$gender,age = df_Mobility_China_na$age,education = df_Mobility_China_na$education,inc_n = df_Mobility_China_na$inc_n,children = df_Mobility_China_na$household_children)),
                                                treat = df_Mobility_China_na$Transport_,
                                                scale.type = "TTX",
                                                baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new3_Mobility_China)

sum3_Mobility_Chinaa<- print(out_new3_Mobility_China)
sum3_Mobility_Chinab <- as.data.frame(sum3_Mobility_Chinaa$statistics)
sum3_Mobility_Chinac <- as.data.frame(sum3_Mobility_Chinaa$quantiles)
sum3_Mobility_China<-cbind (sum3_Mobility_Chinab,sum3_Mobility_Chinac)
sum3_Mobility_China$var <- rownames(sum3_Mobility_China)

sum3_Mobility_China<-sum3_Mobility_China %>%mutate_each(funs(round(., 3)), -var)
sum3_Mobility_China <- sum3_Mobility_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum3_Mobility_China, file = "sum3_Mobility_China.csv", sep = ";", quote = FALSE, row.names = F)

sum3_Mobility_China <- sum3_Mobility_China [1:4,]
sum3_Mobility_China$country <- "China"
sum3_Mobility_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum3_Mobility_China$DV = factor("Willingness to Pay more for Motor Fuels")


sum3_df <- rbind(sum3_Mobility_China,sum3_Mobility_Germany,sum3_Mobility_USA,sum3_food_China,sum3_food_Germany, sum3_food_USA)



#Model 4: DV - Q12.2 Slider/Personal reduction of consumption-------
#Food
out_new4_food_USA <- sparsereg::sparsereg(y = df_food_USA_na$slider_reduce_per_meatcon_z, 
                                          X = as.matrix(data.frame(cri_sust = df_food_USA_na$cri_sust,cri_ego = df_food_USA_na$cri_ego,pv_ego= df_food_USA_na$pv_ego, pv_altr=df_food_USA_na$pv_altr,pv_bios= df_food_USA_na$pv_bios, pv_hedo=df_food_USA_na$pv_hedo,food_habits = df_food_USA_na$food_habits,meals_week = df_food_USA_na$meals_week,amount_meat = df_food_USA_na$amount_meat,dif_reduce = df_food_USA_na$dif_reduce,quality_life = df_food_USA_na$quality_life,occupation = df_food_USA_na$occupation,gender = df_food_USA_na$gender,age = df_food_USA_na$age,education = df_food_USA_na$education,inc_n = df_food_USA_na$inc_n,children = df_food_USA_na$household_children,urbanRural = df_food_USA_na$urbanRural, left_right = df_food_USA_na$left_right, government_int = df_food_USA_na$government_int)),
                                          treat = df_food_USA_na$Food_,
                                          scale.type = "TTX",
                                          baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new4_food_USA)

sum4_food_USAa<- print(out_new4_food_USA)
sum4_food_USAb <- as.data.frame(sum4_food_USAa$statistics)
sum4_food_USAc <- as.data.frame(sum4_food_USAa$quantiles)
sum4_food_USA<-cbind (sum4_food_USAb,sum4_food_USAc)
sum4_food_USA$var <- rownames(sum4_food_USA)

sum4_food_USA<-sum4_food_USA %>%mutate_each(funs(round(., 3)), -var)
sum4_food_USA <- sum4_food_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum4_food_USA, file = "sum4_food_USA.csv", sep = ";", quote = FALSE, row.names = F)


sum4_food_USA <- sum4_food_USA [1:4,]
sum4_food_USA$country <- "USA"
sum4_food_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum4_food_USA$DV = factor("Intentions to Reduce Meat/Fish Consumption")


out_new4_food_Germany <- sparsereg::sparsereg(y = df_food_Germany_na$slider_reduce_per_meatcon_z, 
                                              X = as.matrix(data.frame(cri_sust = df_food_Germany_na$cri_sust,cri_ego = df_food_Germany_na$cri_ego,pv_ego= df_food_Germany_na$pv_ego, pv_altr=df_food_Germany_na$pv_altr,pv_bios= df_food_Germany_na$pv_bios, pv_hedo=df_food_Germany_na$pv_hedo,food_habits = df_food_Germany_na$food_habits,meals_week = df_food_Germany_na$meals_week,amount_meat = df_food_Germany_na$amount_meat,dif_reduce = df_food_Germany_na$dif_reduce,quality_life = df_food_Germany_na$quality_life,occupation = df_food_Germany_na$occupation,gender = df_food_Germany_na$gender,age = df_food_Germany_na$age,education = df_food_Germany_na$education,inc_n = df_food_Germany_na$inc_n,children = df_food_Germany_na$household_children,urbanRural = df_food_Germany_na$urbanRural, left_right = df_food_Germany_na$left_right, government_int = df_food_Germany_na$government_int)),
                                              treat = df_food_Germany_na$Food_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new4_food_Germany)


sum4_food_Germanya<- print(out_new4_food_Germany)
sum4_food_Germanyb <- as.data.frame(sum4_food_Germanya$statistics)
sum4_food_Germanyc <- as.data.frame(sum4_food_Germanya$quantiles)
sum4_food_Germany<-cbind (sum4_food_Germanyb,sum4_food_Germanyc)
sum4_food_Germany$var <- rownames(sum4_food_Germany)

sum4_food_Germany<-sum4_food_Germany %>%mutate_each(funs(round(., 3)), -var)
sum4_food_Germany <- sum4_food_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum4_food_Germany, file = "sum4_food_Germany.csv", sep = ";", quote = FALSE, row.names = F)

sum4_food_Germany <- sum4_food_Germany [1:4,]
sum4_food_Germany$country <- "Germany"
sum4_food_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum4_food_Germany$DV = factor("Intentions to Reduce Meat/Fish Consumption")



out_new4_food_China <- sparsereg::sparsereg(y = df_food_China_na$slider_reduce_per_meatcon_z, 
                                            X = as.matrix(data.frame(cri_sust = df_food_China_na$cri_sust,cri_ego = df_food_China_na$cri_ego,pv_ego= df_food_China_na$pv_ego, pv_altr=df_food_China_na$pv_altr,pv_bios= df_food_China_na$pv_bios, pv_hedo=df_food_China_na$pv_hedo,food_habits = df_food_China_na$food_habits,meals_week = df_food_China_na$meals_week,amount_meat = df_food_China_na$amount_meat,dif_reduce = df_food_China_na$dif_reduce,quality_life = df_food_China_na$quality_life,occupation = df_food_China_na$occupation,gender = df_food_China_na$gender,age = df_food_China_na$age,education = df_food_China_na$education,inc_n = df_food_China_na$inc_n,children = df_food_China_na$household_children)),
                                            treat = df_food_China_na$Food_,
                                            scale.type = "TTX",
                                            baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new4_food_China)


sum4_food_Chinaa<- print(out_new4_food_China)
sum4_food_Chinab <- as.data.frame(sum4_food_Chinaa$statistics)
sum4_food_Chinac <- as.data.frame(sum4_food_Chinaa$quantiles)
sum4_food_China<-cbind (sum4_food_Chinab,sum4_food_Chinac)
sum4_food_China$var <- rownames(sum4_food_China)

sum4_food_China<-sum4_food_China %>%mutate_each(funs(round(., 3)), -var)
sum4_food_China <- sum4_food_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum4_food_China, file = "sum4_food_China.csv", sep = ";", quote = FALSE, row.names = F)

sum4_food_China <- sum4_food_China [1:4,]
sum4_food_China$country <- "China"
sum4_food_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum4_food_China$DV = factor("Intentions to Reduce Meat/Fish Consumption")

#Mobility
out_new4_Mobility_USA <- sparsereg::sparsereg(y = df_Mobility_USA_na$slider_reduce_per_caruse_z, 
                                              X = as.matrix(data.frame(crim_sust = df_Mobility_USA_na$crim_sust,crim_ego = df_Mobility_USA_na$crim_ego,pv_ego= df_Mobility_USA_na$pv_ego, pv_altr=df_Mobility_USA_na$pv_altr,pv_bios= df_Mobility_USA_na$pv_bios, pv_hedo=df_Mobility_USA_na$pv_hedo,driving_habits = df_Mobility_USA_na$driving_habits,drivers_license=df_Mobility_USA_na$drivers_license,car_access= df_Mobility_USA_na$car_access,car_engine=df_Mobility_USA_na$car_engine,  km_driven  = df_Mobility_USA_na$km_driven,dif_drivstop=df_Mobility_USA_na$dif_drivstop, quality_life_m=df_Mobility_USA_na$quality_life_m, occupation = df_Mobility_USA_na$occupation,gender = df_Mobility_USA_na$gender,age = df_Mobility_USA_na$age,education = df_Mobility_USA_na$education,inc_n = df_Mobility_USA_na$inc_n,children = df_Mobility_USA_na$household_children,urbanRural = df_Mobility_USA_na$urbanRural, left_right = df_Mobility_USA_na$left_right, government_int = df_Mobility_USA_na$government_int)),
                                              treat = df_Mobility_USA_na$Transport_,
                                              scale.type = "TTX",
                                              baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new4_Mobility_USA)

sum4_Mobility_USAa<- print(out_new4_Mobility_USA)
sum4_Mobility_USAb <- as.data.frame(sum4_Mobility_USAa$statistics)
sum4_Mobility_USAc <- as.data.frame(sum4_Mobility_USAa$quantiles)
sum4_Mobility_USA<-cbind (sum4_Mobility_USAb,sum4_Mobility_USAc)
sum4_Mobility_USA$var <- rownames(sum4_Mobility_USA)

sum4_Mobility_USA<-sum4_Mobility_USA %>%mutate_each(funs(round(., 3)), -var)
sum4_Mobility_USA <- sum4_Mobility_USA %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum4_Mobility_USA, file = "sum4_Mobility_USA.csv", sep = ";", quote = FALSE, row.names = F)

sum4_Mobility_USA <- sum4_Mobility_USA [1:4,]
sum4_Mobility_USA$country <- "USA"
sum4_Mobility_USA$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum4_Mobility_USA$DV = factor("Intentions to Reduce Car Use")


out_new4_Mobility_Germany <- sparsereg::sparsereg(y = df_Mobility_Germany_na$slider_reduce_per_caruse_z, 
                                                  X = as.matrix(data.frame(crim_sust = df_Mobility_Germany_na$crim_sust,crim_ego = df_Mobility_Germany_na$crim_ego,pv_ego= df_Mobility_Germany_na$pv_ego, pv_altr=df_Mobility_Germany_na$pv_altr,pv_bios= df_Mobility_Germany_na$pv_bios, pv_hedo=df_Mobility_Germany_na$pv_hedo,driving_habits = df_Mobility_Germany_na$driving_habits,drivers_license=df_Mobility_Germany_na$drivers_license,car_access= df_Mobility_Germany_na$car_access,car_engine=df_Mobility_Germany_na$car_engine,  km_driven  = df_Mobility_Germany_na$km_driven,dif_drivstop=df_Mobility_Germany_na$dif_drivstop, quality_life_m=df_Mobility_Germany_na$quality_life_m, occupation = df_Mobility_Germany_na$occupation,gender = df_Mobility_Germany_na$gender,age = df_Mobility_Germany_na$age,education = df_Mobility_Germany_na$education,inc_n = df_Mobility_Germany_na$inc_n,children = df_Mobility_Germany_na$household_children,urbanRural = df_Mobility_Germany_na$urbanRural, left_right = df_Mobility_Germany_na$left_right, government_int = df_Mobility_Germany_na$government_int)),
                                                  treat = df_Mobility_Germany_na$Transport_,
                                                  scale.type = "TTX",
                                                  baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new4_Mobility_Germany)

sum4_Mobility_Germanya<- print(out_new4_Mobility_Germany)
sum4_Mobility_Germanyb <- as.data.frame(sum4_Mobility_Germanya$statistics)
sum4_Mobility_Germanyc <- as.data.frame(sum4_Mobility_Germanya$quantiles)
sum4_Mobility_Germany<-cbind (sum4_Mobility_Germanyb,sum4_Mobility_Germanyc)
sum4_Mobility_Germany$var <- rownames(sum4_Mobility_Germany)

sum4_Mobility_Germany<-sum4_Mobility_Germany %>%mutate_each(funs(round(., 3)), -var)
sum4_Mobility_Germany <- sum4_Mobility_Germany %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum4_Mobility_Germany, file = "sum4_Mobility_Germany.csv", sep = ";", quote = FALSE, row.names = F)

sum4_Mobility_Germany <- sum4_Mobility_Germany [1:4,]
sum4_Mobility_Germany$country <- "Germany"
sum4_Mobility_Germany$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum4_Mobility_Germany$DV = factor("Intentions to Reduce Car Use")

out_new4_Mobility_China <- sparsereg::sparsereg(y = df_Mobility_China_na$slider_reduce_per_caruse_z, 
                                                X = as.matrix(data.frame(crim_sust = df_Mobility_China_na$crim_sust,crim_ego = df_Mobility_China_na$crim_ego,pv_ego= df_Mobility_China_na$pv_ego, pv_altr=df_Mobility_China_na$pv_altr,pv_bios= df_Mobility_China_na$pv_bios, pv_hedo=df_Mobility_China_na$pv_hedo,driving_habits = df_Mobility_China_na$driving_habits,drivers_license=df_Mobility_China_na$drivers_license,car_access= df_Mobility_China_na$car_access,car_engine=df_Mobility_China_na$car_engine,  km_driven  = df_Mobility_China_na$km_driven,dif_drivstop=df_Mobility_China_na$dif_drivstop, quality_life_m=df_Mobility_China_na$quality_life_m, occupation = df_Mobility_China_na$occupation,gender = df_Mobility_China_na$gender,age = df_Mobility_China_na$age,education = df_Mobility_China_na$education,inc_n = df_Mobility_China_na$inc_n,children = df_Mobility_China_na$household_children)),
                                                treat = df_Mobility_China_na$Transport_,
                                                scale.type = "TTX",
                                                baseline.vec = c("Control"))

sparsereg::summary.sparsereg(out_new4_Mobility_China)

sum4_Mobility_Chinaa<- print(out_new4_Mobility_China)
sum4_Mobility_Chinab <- as.data.frame(sum4_Mobility_Chinaa$statistics)
sum4_Mobility_Chinac <- as.data.frame(sum4_Mobility_Chinaa$quantiles)
sum4_Mobility_China<-cbind (sum4_Mobility_Chinab,sum4_Mobility_Chinac)
sum4_Mobility_China$var <- rownames(sum4_Mobility_China)

sum4_Mobility_China<-sum4_Mobility_China %>%mutate_each(funs(round(., 3)), -var)
sum4_Mobility_China <- sum4_Mobility_China %>% select(var, Mean,`50%`,`2.5%`,`97.5%`)
write.table(sum4_Mobility_China, file = "sum4_Mobility_China.csv", sep = ";", quote = FALSE, row.names = F)

sum4_Mobility_China <- sum4_Mobility_China [1:4,]
sum4_Mobility_China$country <- "China"
sum4_Mobility_China$Treatment <- factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth"))
sum4_Mobility_China$DV = factor("Intentions to Reduce Car Use")


sum4_df <- rbind(sum4_Mobility_China,sum4_Mobility_Germany,sum4_Mobility_USA,sum4_food_China,sum4_food_Germany, sum4_food_USA)


#Estimation of credible intervals around the median, see details https://www.sciencedirect.com/topics/mathematics/credible-interval

#Linear Regression Models Restricted Sample (Only individuals that eat meat/fish or drive a car themselves)---------
#Model 1 - DV - Concern about Impact of Meat/Fish Consumption/Car Use----------

#Food
M1_food_usa_ITT <- lm_robust(concern_z ~ Food_ , data= df_food_USA_na )
summary(M1_food_usa_ITT)

M1_food_germany_ITT <- lm_robust(concern_z ~ Food_ , data= df_food_Germany_na )
summary(M1_food_germany_ITT)

M1_food_china_ITT <- lm_robust(concern_z ~ Food_ , data= df_food_China_na)
summary(M1_food_china_ITT)

#Preparing Plots for Main Effects
mydf_M1_food_usa_ITT <- data.frame(
  coef = c(coef(summary(M1_food_usa_ITT ))[2],coef(summary(M1_food_usa_ITT))[3], coef(summary(M1_food_usa_ITT))[4], coef(summary(M1_food_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_food_usa_ITT ))[2,2], coef(summary(M1_food_usa_ITT))[3,2], coef(summary(M1_food_usa_ITT))[4,2], coef(summary(M1_food_usa_ITT))[5,2]),
  DV = factor(rep(c("Concern about Meat/Fish Consumption"), 4)),
  country = "USA")

mydf_M1_food_germany_ITT <- data.frame(
  coef = c(coef(summary(M1_food_germany_ITT))[2],coef(summary(M1_food_germany_ITT))[3], coef(summary(M1_food_germany_ITT))[4], coef(summary(M1_food_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_food_germany_ITT ))[2,2], coef(summary(M1_food_germany_ITT))[3,2], coef(summary(M1_food_germany_ITT))[4,2], coef(summary(M1_food_germany_ITT))[5,2]),
  DV = factor(rep(c("Concern about Meat/Fish Consumption"), 4)),
  country = "Germany")

mydf_M1_food_china_ITT <- data.frame(
  coef = c(coef(summary(M1_food_china_ITT))[2],coef(summary(M1_food_china_ITT))[3], coef(summary(M1_food_china_ITT))[4], coef(summary(M1_food_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_food_china_ITT ))[2,2], coef(summary(M1_food_china_ITT))[3,2], coef(summary(M1_food_china_ITT))[4,2], coef(summary(M1_food_china_ITT))[5,2]),
  DV = factor(rep(c("Concern about Meat/Fish Consumption"), 4)),
  country = "China")


b1 <- rbind(mydf_M1_food_china_ITT,mydf_M1_food_germany_ITT, mydf_M1_food_usa_ITT)

b1 <- b1[, names(b1) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

#Mobility

M1_mobility_usa_ITT <- lm_robust(concern_z ~ treatment_mobility , data=df_Mobility_USA_na)
summary(M1_mobility_usa_ITT)

M1_mobility_germany_ITT <- lm_robust(concern_z ~ treatment_mobility , data= df_Mobility_Germany_na)
summary(M1_mobility_germany_ITT)

M1_mobility_china_ITT <- lm_robust(concern_z ~ treatment_mobility , data= df_Mobility_China_na)
summary(M1_mobility_china_ITT)

#Preparing Plots for Main Effects

mydf_M1_mobility_usa_ITT <- data.frame(
  coef = c(coef(summary(M1_mobility_usa_ITT))[2],coef(summary(M1_mobility_usa_ITT))[3], coef(summary(M1_mobility_usa_ITT))[4], coef(summary(M1_mobility_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_mobility_usa_ITT ))[2,2], coef(summary(M1_mobility_usa_ITT))[3,2], coef(summary(M1_mobility_usa_ITT))[4,2], coef(summary(M1_mobility_usa_ITT))[5,2]),
  DV = factor(rep(c("Concern about Car Use"), 4)),
  country = "USA")

mydf_M1_mobility_germany_ITT <- data.frame(
  coef = c(coef(summary(M1_mobility_germany_ITT))[2],coef(summary(M1_mobility_germany_ITT))[3], coef(summary(M1_mobility_germany_ITT))[4], coef(summary(M1_mobility_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_mobility_germany_ITT ))[2,2], coef(summary(M1_mobility_germany_ITT))[3,2], coef(summary(M1_mobility_germany_ITT))[4,2], coef(summary(M1_mobility_germany_ITT))[5,2]),
  DV = factor(rep(c("Concern about Car Use"), 4)),
  country = "Germany")

mydf_M1_mobility_china_ITT <- data.frame(
  coef = c(coef(summary(M1_mobility_china_ITT))[2],coef(summary(M1_mobility_china_ITT))[3], coef(summary(M1_mobility_china_ITT))[4], coef(summary(M1_mobility_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_mobility_china_ITT ))[2,2], coef(summary(M1_mobility_china_ITT))[3,2], coef(summary(M1_mobility_china_ITT))[4,2], coef(summary(M1_mobility_china_ITT))[5,2]),
  DV = factor(rep(c("Concern about Car Use"), 4)),
  country = "China")


c1 <- rbind(mydf_M1_mobility_china_ITT,mydf_M1_mobility_germany_ITT, mydf_M1_mobility_usa_ITT)
c1 <- c1[, names(c1) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d1 <- rbind (c1,b1)

#Model 2: DV - Q14.3.	SA PER ITEM/ Policy Support Reduction----------

M2_food_usa_ITT <- lm_robust(policysupport_reduction_meatcon_z ~ Food_ , data= df_food_USA_na)
summary(M2_food_usa_ITT)

M2_food_germany_ITT <- lm_robust(policysupport_reduction_meatcon_z ~ Food_ , data= df_food_Germany_na)
summary(M2_food_germany_ITT)

M2_food_china_ITT <- lm_robust(policysupport_reduction_meatcon_z ~ Food_ , data= df_food_China_na)
summary(M2_food_china_ITT)

#Preparing Plots for Main Effects
mydf_M2_food_usa_ITT <- data.frame(
  coef = c(coef(summary(M2_food_usa_ITT ))[2],coef(summary(M2_food_usa_ITT))[3], coef(summary(M2_food_usa_ITT))[4], coef(summary(M2_food_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_food_usa_ITT ))[2,2], coef(summary(M2_food_usa_ITT))[3,2], coef(summary(M2_food_usa_ITT))[4,2], coef(summary(M2_food_usa_ITT))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Meat/Fish Consumption"), 4)),
  country = "USA")

mydf_M2_food_germany_ITT <- data.frame(
  coef = c(coef(summary(M2_food_germany_ITT))[2],coef(summary(M2_food_germany_ITT))[3], coef(summary(M2_food_germany_ITT))[4], coef(summary(M2_food_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_food_germany_ITT ))[2,2], coef(summary(M2_food_germany_ITT))[3,2], coef(summary(M2_food_germany_ITT))[4,2], coef(summary(M2_food_germany_ITT))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Meat/Fish Consumption"), 4)),
  country = "Germany")

mydf_M2_food_china_ITT <- data.frame(
  coef = c(coef(summary(M2_food_china_ITT))[2],coef(summary(M2_food_china_ITT))[3], coef(summary(M2_food_china_ITT))[4], coef(summary(M2_food_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_food_china_ITT ))[2,2], coef(summary(M2_food_china_ITT))[3,2], coef(summary(M2_food_china_ITT))[4,2], coef(summary(M2_food_china_ITT))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Meat/Fish Consumption"), 4)),
  country = "China")


b2 <- rbind(mydf_M2_food_china_ITT,mydf_M2_food_germany_ITT, mydf_M2_food_usa_ITT)

b2 <- b2[, names(b2) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

#Mobility

M2_mobility_usa_ITT <- lm_robust(policysupport_reduction_caruse_z ~ treatment_mobility , data= df_Mobility_USA_na)
summary(M2_mobility_usa_ITT)

M2_mobility_germany_ITT <- lm_robust(policysupport_reduction_caruse_z ~ treatment_mobility , data= df_Mobility_Germany_na)
summary(M2_mobility_germany_ITT)

M2_mobility_china_ITT <- lm_robust(policysupport_reduction_caruse_z ~ treatment_mobility , data= df_Mobility_China_na)
summary(M2_mobility_china_ITT)

#Preparing Plots for Main Effects

mydf_M2_mobility_usa_ITT <- data.frame(
  coef = c(coef(summary(M2_mobility_usa_ITT))[2],coef(summary(M2_mobility_usa_ITT))[3], coef(summary(M2_mobility_usa_ITT))[4], coef(summary(M2_mobility_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_mobility_usa_ITT ))[2,2], coef(summary(M2_mobility_usa_ITT))[3,2], coef(summary(M2_mobility_usa_ITT))[4,2], coef(summary(M2_mobility_usa_ITT))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Car Use"), 4)),
  country = "USA")

mydf_M2_mobility_germany_ITT <- data.frame(
  coef = c(coef(summary(M2_mobility_germany_ITT))[2],coef(summary(M2_mobility_germany_ITT))[3], coef(summary(M2_mobility_germany_ITT))[4], coef(summary(M2_mobility_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_mobility_germany_ITT ))[2,2], coef(summary(M2_mobility_germany_ITT))[3,2], coef(summary(M2_mobility_germany_ITT))[4,2], coef(summary(M2_mobility_germany_ITT))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Car Use"), 4)),
  country = "Germany")

mydf_M2_mobility_china_ITT <- data.frame(
  coef = c(coef(summary(M2_mobility_china_ITT))[2],coef(summary(M2_mobility_china_ITT))[3], coef(summary(M2_mobility_china_ITT))[4], coef(summary(M2_mobility_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_mobility_china_ITT ))[2,2], coef(summary(M2_mobility_china_ITT))[3,2], coef(summary(M2_mobility_china_ITT))[4,2], coef(summary(M2_mobility_china_ITT))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Car Use"), 4)),
  country = "China")


c2 <- rbind(mydf_M2_mobility_china_ITT,mydf_M2_mobility_germany_ITT, mydf_M2_mobility_usa_ITT)
c2 <- c2[, names(c2) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d2 <- rbind (c2,b2)

#Model 3: DV - Q14.1 Slider/Personal WTP--------
M3_food_usa_ITT <- lm_robust(slider_WTP_meat_z ~ Food_, data= df_food_USA_na)
summary(M3_food_usa_ITT)

M3_food_germany_ITT <- lm_robust(slider_WTP_meat_z ~ Food_ , data= df_food_Germany_na)
summary(M3_food_germany_ITT)

M3_food_china_ITT <- lm_robust(slider_WTP_meat_z ~ Food_ , data= df_food_China_na)
summary(M3_food_china_ITT)

#Preparing Plots for Main Effects
mydf_M3_food_usa_ITT <- data.frame(
  coef = c(coef(summary(M3_food_usa_ITT ))[2],coef(summary(M3_food_usa_ITT))[3], coef(summary(M3_food_usa_ITT))[4], coef(summary(M3_food_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_food_usa_ITT ))[2,2], coef(summary(M3_food_usa_ITT))[3,2], coef(summary(M3_food_usa_ITT))[4,2], coef(summary(M3_food_usa_ITT))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Meat/Fish"), 4)),
  country = "USA")

mydf_M3_food_germany_ITT <- data.frame(
  coef = c(coef(summary(M3_food_germany_ITT))[2],coef(summary(M3_food_germany_ITT))[3], coef(summary(M3_food_germany_ITT))[4], coef(summary(M3_food_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_food_germany_ITT ))[2,2], coef(summary(M3_food_germany_ITT))[3,2], coef(summary(M3_food_germany_ITT))[4,2], coef(summary(M3_food_germany_ITT))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Meat/Fish"), 4)),
  country = "Germany")

mydf_M3_food_china_ITT <- data.frame(
  coef = c(coef(summary(M3_food_china_ITT))[2],coef(summary(M3_food_china_ITT))[3], coef(summary(M3_food_china_ITT))[4], coef(summary(M3_food_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_food_china_ITT ))[2,2], coef(summary(M3_food_china_ITT))[3,2], coef(summary(M3_food_china_ITT))[4,2], coef(summary(M3_food_china_ITT))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Meat/Fish"), 4)),
  country = "China")


b3 <- rbind(mydf_M3_food_china_ITT,mydf_M3_food_germany_ITT, mydf_M3_food_usa_ITT)

b3 <- b3[, names(b3) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

#Mobility

M3_mobility_usa_ITT <- lm_robust(slider_WTP_mobility_z ~ treatment_mobility , data= df_Mobility_USA_na)
summary(M3_mobility_usa_ITT)

M3_mobility_germany_ITT <- lm_robust(slider_WTP_mobility_z ~ treatment_mobility , data= df_Mobility_Germany_na)
summary(M3_mobility_germany_ITT)

M3_mobility_china_ITT <- lm_robust(slider_WTP_mobility_z ~ treatment_mobility , data= df_Mobility_China_na)
summary(M3_mobility_china_ITT)

#Preparing Plots for Main Effects

mydf_M3_mobility_usa_ITT <- data.frame(
  coef = c(coef(summary(M3_mobility_usa_ITT))[2],coef(summary(M3_mobility_usa_ITT))[3], coef(summary(M3_mobility_usa_ITT))[4], coef(summary(M3_mobility_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_mobility_usa_ITT ))[2,2], coef(summary(M3_mobility_usa_ITT))[3,2], coef(summary(M3_mobility_usa_ITT))[4,2], coef(summary(M3_mobility_usa_ITT))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Motor Fuels"), 4)),
  country = "USA")

mydf_M3_mobility_germany_ITT <- data.frame(
  coef = c(coef(summary(M3_mobility_germany_ITT))[2],coef(summary(M3_mobility_germany_ITT))[3], coef(summary(M3_mobility_germany_ITT))[4], coef(summary(M3_mobility_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_mobility_germany_ITT ))[2,2], coef(summary(M3_mobility_germany_ITT))[3,2], coef(summary(M3_mobility_germany_ITT))[4,2], coef(summary(M3_mobility_germany_ITT))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Motor Fuels"), 4)),
  country = "Germany")

mydf_M3_mobility_china_ITT <- data.frame(
  coef = c(coef(summary(M3_mobility_china_ITT))[2],coef(summary(M3_mobility_china_ITT))[3], coef(summary(M3_mobility_china_ITT))[4], coef(summary(M3_mobility_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_mobility_china_ITT ))[2,2], coef(summary(M3_mobility_china_ITT))[3,2], coef(summary(M3_mobility_china_ITT))[4,2], coef(summary(M3_mobility_china_ITT))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Motor Fuels"), 4)),
  country = "China")


c3 <- rbind(mydf_M3_mobility_china_ITT,mydf_M3_mobility_germany_ITT, mydf_M3_mobility_usa_ITT)
c3 <- c3[, names(c3) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d3<- rbind (c3,b3)


#Model 4: DV - Q12.2 Slider/Personal reduction of consumption-------

#Food
M4_food_usa_ITT <- lm_robust(slider_reduce_per_meatcon_z ~ Food_ , data= df_food_USA_na)
summary(M4_food_usa_ITT)

M4_food_germany_ITT <- lm_robust(slider_reduce_per_meatcon_z ~ Food_ , data= df_food_Germany_na)
summary(M4_food_germany_ITT)

M4_food_china_ITT <- lm_robust(slider_reduce_per_meatcon_z ~ Food_ , data= df_food_China_na)
summary(M4_food_china_ITT)

#Preparing Plots for Main Effects
mydf_M4_food_usa_ITT <- data.frame(
  coef = c(coef(summary(M4_food_usa_ITT ))[2],coef(summary(M4_food_usa_ITT))[3], coef(summary(M4_food_usa_ITT))[4], coef(summary(M4_food_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_food_usa_ITT ))[2,2], coef(summary(M4_food_usa_ITT))[3,2], coef(summary(M4_food_usa_ITT))[4,2], coef(summary(M4_food_usa_ITT))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Meat/Fish Consumption"), 4)),
  country = "USA")

mydf_M4_food_germany_ITT <- data.frame(
  coef = c(coef(summary(M4_food_germany_ITT))[2],coef(summary(M4_food_germany_ITT))[3], coef(summary(M4_food_germany_ITT))[4], coef(summary(M4_food_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_food_germany_ITT ))[2,2], coef(summary(M4_food_germany_ITT))[3,2], coef(summary(M4_food_germany_ITT))[4,2], coef(summary(M4_food_germany_ITT))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Meat/Fish Consumption"), 4)),
  country = "Germany")

mydf_M4_food_china_ITT <- data.frame(
  coef = c(coef(summary(M4_food_china_ITT))[2],coef(summary(M4_food_china_ITT))[3], coef(summary(M4_food_china_ITT))[4], coef(summary(M4_food_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_food_china_ITT ))[2,2], coef(summary(M4_food_china_ITT))[3,2], coef(summary(M4_food_china_ITT))[4,2], coef(summary(M4_food_china_ITT))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Meat/Fish Consumption"), 4)),
  country = "China")


b4 <- rbind(mydf_M4_food_china_ITT,mydf_M4_food_germany_ITT, mydf_M4_food_usa_ITT)

b4 <- b4[, names(b4) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 


#Mobility

M4_mobility_usa_ITT <- lm_robust(slider_reduce_per_caruse_z ~ treatment_mobility , data= df_Mobility_USA_na)
summary(M4_mobility_usa_ITT)

M4_mobility_germany_ITT <- lm_robust(slider_reduce_per_caruse_z ~ treatment_mobility , data= df_Mobility_Germany_na)
summary(M4_mobility_germany_ITT)

M4_mobility_china_ITT <- lm_robust(slider_reduce_per_caruse_z ~ treatment_mobility , data= df_Mobility_China_na)
summary(M4_mobility_china_ITT)

#Preparing Plots for Main Effects

mydf_M4_mobility_usa_ITT <- data.frame(
  coef = c(coef(summary(M4_mobility_usa_ITT))[2],coef(summary(M4_mobility_usa_ITT))[3], coef(summary(M4_mobility_usa_ITT))[4], coef(summary(M4_mobility_usa_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_mobility_usa_ITT ))[2,2], coef(summary(M4_mobility_usa_ITT))[3,2], coef(summary(M4_mobility_usa_ITT))[4,2], coef(summary(M4_mobility_usa_ITT))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Car Use"), 4)),
  country = "USA")

mydf_M4_mobility_germany_ITT <- data.frame(
  coef = c(coef(summary(M4_mobility_germany_ITT))[2],coef(summary(M4_mobility_germany_ITT))[3], coef(summary(M4_mobility_germany_ITT))[4], coef(summary(M4_mobility_germany_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_mobility_germany_ITT ))[2,2], coef(summary(M4_mobility_germany_ITT))[3,2], coef(summary(M4_mobility_germany_ITT))[4,2], coef(summary(M4_mobility_germany_ITT))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Car Use"), 4)),
  country = "Germany")

mydf_M4_mobility_china_ITT <- data.frame(
  coef = c(coef(summary(M4_mobility_china_ITT))[2],coef(summary(M4_mobility_china_ITT))[3], coef(summary(M4_mobility_china_ITT))[4], coef(summary(M4_mobility_china_ITT ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_mobility_china_ITT ))[2,2], coef(summary(M4_mobility_china_ITT))[3,2], coef(summary(M4_mobility_china_ITT))[4,2], coef(summary(M4_mobility_china_ITT))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Car Use"), 4)),
  country = "China")


c4 <- rbind(mydf_M4_mobility_china_ITT,mydf_M4_mobility_germany_ITT, mydf_M4_mobility_usa_ITT)
c4<- c4[, names(c4) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d4 <- rbind (c4,b4)


#Robustness Check: Linear Regression Models Full Sample---------
#Model 1 - DV - Concern about Impact of Meat/Fish Consumption/Car Use----------

#Food
M1_food_usa_ITT_Robust <- lm_robust(concern_z ~ Food_ , data= subset (df, group =="Food"& country=="USA" ))
summary(M1_food_usa_ITT_Robust)

M1_food_germany_ITT_Robust <- lm_robust(concern_z ~ Food_ , data= subset (df, group =="Food"& country=="Germany" ))
summary(M1_food_germany_ITT_Robust)

M1_food_china_ITT_Robust <- lm_robust(concern_z ~ Food_ , data= subset (df, group =="Food"& country=="China" ))
summary(M1_food_china_ITT_Robust)

#Preparing Plots for Main Effects
mydf_M1_food_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M1_food_usa_ITT_Robust ))[2],coef(summary(M1_food_usa_ITT_Robust))[3], coef(summary(M1_food_usa_ITT_Robust))[4], coef(summary(M1_food_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_food_usa_ITT_Robust ))[2,2], coef(summary(M1_food_usa_ITT_Robust))[3,2], coef(summary(M1_food_usa_ITT_Robust))[4,2], coef(summary(M1_food_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Concern about Meat/Fish Consumption"), 4)),
  country = "USA")

mydf_M1_food_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M1_food_germany_ITT_Robust))[2],coef(summary(M1_food_germany_ITT_Robust))[3], coef(summary(M1_food_germany_ITT_Robust))[4], coef(summary(M1_food_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_food_germany_ITT_Robust ))[2,2], coef(summary(M1_food_germany_ITT_Robust))[3,2], coef(summary(M1_food_germany_ITT_Robust))[4,2], coef(summary(M1_food_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Concern about Meat/Fish Consumption"), 4)),
  country = "Germany")

mydf_M1_food_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M1_food_china_ITT_Robust))[2],coef(summary(M1_food_china_ITT_Robust))[3], coef(summary(M1_food_china_ITT_Robust))[4], coef(summary(M1_food_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_food_china_ITT_Robust ))[2,2], coef(summary(M1_food_china_ITT_Robust))[3,2], coef(summary(M1_food_china_ITT_Robust))[4,2], coef(summary(M1_food_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Concern about Meat/Fish Consumption"), 4)),
  country = "China")


b1_Robust <- rbind(mydf_M1_food_china_ITT_Robust,mydf_M1_food_germany_ITT_Robust, mydf_M1_food_usa_ITT_Robust)

b1_Robust <- b1_Robust[, names(b1_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

#Mobility

M1_mobility_usa_ITT_Robust <- lm_robust(concern_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="USA"))
summary(M1_mobility_usa_ITT_Robust)

M1_mobility_germany_ITT_Robust <- lm_robust(concern_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="Germany"))
summary(M1_mobility_germany_ITT_Robust)

M1_mobility_china_ITT_Robust <- lm_robust(concern_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="China"))
summary(M1_mobility_china_ITT_Robust)

#Preparing Plots for Main Effects

mydf_M1_mobility_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M1_mobility_usa_ITT_Robust))[2],coef(summary(M1_mobility_usa_ITT_Robust))[3], coef(summary(M1_mobility_usa_ITT_Robust))[4], coef(summary(M1_mobility_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_mobility_usa_ITT_Robust ))[2,2], coef(summary(M1_mobility_usa_ITT_Robust))[3,2], coef(summary(M1_mobility_usa_ITT_Robust))[4,2], coef(summary(M1_mobility_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Concern about Car Use"), 4)),
  country = "USA")

mydf_M1_mobility_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M1_mobility_germany_ITT_Robust))[2],coef(summary(M1_mobility_germany_ITT_Robust))[3], coef(summary(M1_mobility_germany_ITT_Robust))[4], coef(summary(M1_mobility_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_mobility_germany_ITT_Robust ))[2,2], coef(summary(M1_mobility_germany_ITT_Robust))[3,2], coef(summary(M1_mobility_germany_ITT_Robust))[4,2], coef(summary(M1_mobility_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Concern about Car Use"), 4)),
  country = "Germany")

mydf_M1_mobility_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M1_mobility_china_ITT_Robust))[2],coef(summary(M1_mobility_china_ITT_Robust))[3], coef(summary(M1_mobility_china_ITT_Robust))[4], coef(summary(M1_mobility_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M1_mobility_china_ITT_Robust ))[2,2], coef(summary(M1_mobility_china_ITT_Robust))[3,2], coef(summary(M1_mobility_china_ITT_Robust))[4,2], coef(summary(M1_mobility_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Concern about Car Use"), 4)),
  country = "China")


c1_Robust <- rbind(mydf_M1_mobility_china_ITT_Robust,mydf_M1_mobility_germany_ITT_Robust, mydf_M1_mobility_usa_ITT_Robust)
c1_Robust <- c1_Robust[, names(c1_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d1_Robust <- rbind (c1_Robust,b1_Robust)

#Model 2: DV - Q14.3.	SA PER ITEM/ Policy Support Reduction----------

M2_food_usa_ITT_Robust <- lm_robust(policysupport_reduction_meatcon_z ~ Food_ , data= subset (df,group =="Food"& country=="USA"))
summary(M2_food_usa_ITT_Robust)

M2_food_germany_ITT_Robust <- lm_robust(policysupport_reduction_meatcon_z ~ Food_ , data= subset (df, group =="Food"& country=="Germany"))
summary(M2_food_germany_ITT_Robust)

M2_food_china_ITT_Robust <- lm_robust(policysupport_reduction_meatcon_z ~ Food_ , data= subset (df, group =="Food"& country=="China"))
summary(M2_food_china_ITT_Robust)

#Preparing Plots for Main Effects
mydf_M2_food_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M2_food_usa_ITT_Robust ))[2],coef(summary(M2_food_usa_ITT_Robust))[3], coef(summary(M2_food_usa_ITT_Robust))[4], coef(summary(M2_food_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_food_usa_ITT_Robust ))[2,2], coef(summary(M2_food_usa_ITT_Robust))[3,2], coef(summary(M2_food_usa_ITT_Robust))[4,2], coef(summary(M2_food_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Meat/Fish Consumption"), 4)),
  country = "USA")

mydf_M2_food_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M2_food_germany_ITT_Robust))[2],coef(summary(M2_food_germany_ITT_Robust))[3], coef(summary(M2_food_germany_ITT_Robust))[4], coef(summary(M2_food_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_food_germany_ITT_Robust ))[2,2], coef(summary(M2_food_germany_ITT_Robust))[3,2], coef(summary(M2_food_germany_ITT_Robust))[4,2], coef(summary(M2_food_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Meat/Fish Consumption"), 4)),
  country = "Germany")

mydf_M2_food_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M2_food_china_ITT_Robust))[2],coef(summary(M2_food_china_ITT_Robust))[3], coef(summary(M2_food_china_ITT_Robust))[4], coef(summary(M2_food_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_food_china_ITT_Robust ))[2,2], coef(summary(M2_food_china_ITT_Robust))[3,2], coef(summary(M2_food_china_ITT_Robust))[4,2], coef(summary(M2_food_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Meat/Fish Consumption"), 4)),
  country = "China")


b2_Robust <- rbind(mydf_M2_food_china_ITT_Robust,mydf_M2_food_germany_ITT_Robust, mydf_M2_food_usa_ITT_Robust)

b2_Robust <- b2_Robust[, names(b2_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

#Mobility

M2_mobility_usa_ITT_Robust <- lm_robust(policysupport_reduction_caruse_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="USA"))
summary(M2_mobility_usa_ITT_Robust)

M2_mobility_germany_ITT_Robust <- lm_robust(policysupport_reduction_caruse_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="Germany" ))
summary(M2_mobility_germany_ITT_Robust)

M2_mobility_china_ITT_Robust <- lm_robust(policysupport_reduction_caruse_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="China"))
summary(M2_mobility_china_ITT_Robust)

#Preparing Plots for Main Effects

mydf_M2_mobility_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M2_mobility_usa_ITT_Robust))[2],coef(summary(M2_mobility_usa_ITT_Robust))[3], coef(summary(M2_mobility_usa_ITT_Robust))[4], coef(summary(M2_mobility_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_mobility_usa_ITT_Robust ))[2,2], coef(summary(M2_mobility_usa_ITT_Robust))[3,2], coef(summary(M2_mobility_usa_ITT_Robust))[4,2], coef(summary(M2_mobility_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Car Use"), 4)),
  country = "USA")

mydf_M2_mobility_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M2_mobility_germany_ITT_Robust))[2],coef(summary(M2_mobility_germany_ITT_Robust))[3], coef(summary(M2_mobility_germany_ITT_Robust))[4], coef(summary(M2_mobility_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_mobility_germany_ITT_Robust ))[2,2], coef(summary(M2_mobility_germany_ITT_Robust))[3,2], coef(summary(M2_mobility_germany_ITT_Robust))[4,2], coef(summary(M2_mobility_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Car Use"), 4)),
  country = "Germany")

mydf_M2_mobility_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M2_mobility_china_ITT_Robust))[2],coef(summary(M2_mobility_china_ITT_Robust))[3], coef(summary(M2_mobility_china_ITT_Robust))[4], coef(summary(M2_mobility_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M2_mobility_china_ITT_Robust ))[2,2], coef(summary(M2_mobility_china_ITT_Robust))[3,2], coef(summary(M2_mobility_china_ITT_Robust))[4,2], coef(summary(M2_mobility_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Policy Support to Reduce Car Use"), 4)),
  country = "China")


c2_Robust <- rbind(mydf_M2_mobility_china_ITT_Robust,mydf_M2_mobility_germany_ITT_Robust, mydf_M2_mobility_usa_ITT_Robust)
c2_Robust <- c2_Robust[, names(c2_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d2_Robust <- rbind (c2_Robust,b2_Robust)

#Model 3: DV - Q14.1 Slider/Personal WTP--------
M3_food_usa_ITT_Robust <- lm_robust(slider_WTP_meat_z ~ Food_, data= subset (df, group =="Food"& country=="USA" ))
summary(M3_food_usa_ITT_Robust)

M3_food_germany_ITT_Robust <- lm_robust(slider_WTP_meat_z ~ Food_ , data= subset (df, group =="Food"& country=="Germany" ))
summary(M3_food_germany_ITT_Robust)

M3_food_china_ITT_Robust <- lm_robust(slider_WTP_meat_z ~ Food_ , data= subset (df, group =="Food"& country=="China" ))
summary(M3_food_china_ITT_Robust)

#Preparing Plots for Main Effects
mydf_M3_food_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M3_food_usa_ITT_Robust ))[2],coef(summary(M3_food_usa_ITT_Robust))[3], coef(summary(M3_food_usa_ITT_Robust))[4], coef(summary(M3_food_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_food_usa_ITT_Robust ))[2,2], coef(summary(M3_food_usa_ITT_Robust))[3,2], coef(summary(M3_food_usa_ITT_Robust))[4,2], coef(summary(M3_food_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Meat/Fish"), 4)),
  country = "USA")

mydf_M3_food_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M3_food_germany_ITT_Robust))[2],coef(summary(M3_food_germany_ITT_Robust))[3], coef(summary(M3_food_germany_ITT_Robust))[4], coef(summary(M3_food_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_food_germany_ITT_Robust ))[2,2], coef(summary(M3_food_germany_ITT_Robust))[3,2], coef(summary(M3_food_germany_ITT_Robust))[4,2], coef(summary(M3_food_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Meat/Fish"), 4)),
  country = "Germany")

mydf_M3_food_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M3_food_china_ITT_Robust))[2],coef(summary(M3_food_china_ITT_Robust))[3], coef(summary(M3_food_china_ITT_Robust))[4], coef(summary(M3_food_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_food_china_ITT_Robust ))[2,2], coef(summary(M3_food_china_ITT_Robust))[3,2], coef(summary(M3_food_china_ITT_Robust))[4,2], coef(summary(M3_food_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Meat/Fish"), 4)),
  country = "China")


b3_Robust <- rbind(mydf_M3_food_china_ITT_Robust,mydf_M3_food_germany_ITT_Robust, mydf_M3_food_usa_ITT_Robust)

b3_Robust <- b3_Robust[, names(b3_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

#Mobility

M3_mobility_usa_ITT_Robust <- lm_robust(slider_WTP_mobility_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="USA"))
summary(M3_mobility_usa_ITT_Robust)

M3_mobility_germany_ITT_Robust <- lm_robust(slider_WTP_mobility_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="Germany" ))
summary(M3_mobility_germany_ITT_Robust)

M3_mobility_china_ITT_Robust <- lm_robust(slider_WTP_mobility_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="China"))
summary(M3_mobility_china_ITT_Robust)

#Preparing Plots for Main Effects

mydf_M3_mobility_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M3_mobility_usa_ITT_Robust))[2],coef(summary(M3_mobility_usa_ITT_Robust))[3], coef(summary(M3_mobility_usa_ITT_Robust))[4], coef(summary(M3_mobility_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_mobility_usa_ITT_Robust ))[2,2], coef(summary(M3_mobility_usa_ITT_Robust))[3,2], coef(summary(M3_mobility_usa_ITT_Robust))[4,2], coef(summary(M3_mobility_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Motor Fuels"), 4)),
  country = "USA")

mydf_M3_mobility_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M3_mobility_germany_ITT_Robust))[2],coef(summary(M3_mobility_germany_ITT_Robust))[3], coef(summary(M3_mobility_germany_ITT_Robust))[4], coef(summary(M3_mobility_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_mobility_germany_ITT_Robust ))[2,2], coef(summary(M3_mobility_germany_ITT_Robust))[3,2], coef(summary(M3_mobility_germany_ITT_Robust))[4,2], coef(summary(M3_mobility_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Motor Fuels"), 4)),
  country = "Germany")

mydf_M3_mobility_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M3_mobility_china_ITT_Robust))[2],coef(summary(M3_mobility_china_ITT_Robust))[3], coef(summary(M3_mobility_china_ITT_Robust))[4], coef(summary(M3_mobility_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M3_mobility_china_ITT_Robust ))[2,2], coef(summary(M3_mobility_china_ITT_Robust))[3,2], coef(summary(M3_mobility_china_ITT_Robust))[4,2], coef(summary(M3_mobility_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Willingness to Pay more for Motor Fuels"), 4)),
  country = "China")


c3_Robust <- rbind(mydf_M3_mobility_china_ITT_Robust,mydf_M3_mobility_germany_ITT_Robust, mydf_M3_mobility_usa_ITT_Robust)
c3_Robust <- c3_Robust[, names(c3_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d3_Robust<- rbind (c3_Robust,b3_Robust)


#Model 4: DV - Q12.2 Slider/Personal reduction of consumption-------

#Food
M4_food_usa_ITT_Robust <- lm_robust(slider_reduce_per_meatcon_z ~ Food_ , data= subset (df, group =="Food"& country=="USA"  ))
summary(M4_food_usa_ITT_Robust)

M4_food_germany_ITT_Robust <- lm_robust(slider_reduce_per_meatcon_z ~ Food_ , data= subset (df, group =="Food"& country=="Germany" ))
summary(M4_food_germany_ITT_Robust)

M4_food_china_ITT_Robust <- lm_robust(slider_reduce_per_meatcon_z ~ Food_ , data= subset (df, group =="Food"& country=="China"  ))
summary(M4_food_china_ITT_Robust)

#Preparing Plots for Main Effects
mydf_M4_food_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M4_food_usa_ITT_Robust ))[2],coef(summary(M4_food_usa_ITT_Robust))[3], coef(summary(M4_food_usa_ITT_Robust))[4], coef(summary(M4_food_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_food_usa_ITT_Robust ))[2,2], coef(summary(M4_food_usa_ITT_Robust))[3,2], coef(summary(M4_food_usa_ITT_Robust))[4,2], coef(summary(M4_food_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Meat/Fish Consumption"), 4)),
  country = "USA")

mydf_M4_food_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M4_food_germany_ITT_Robust))[2],coef(summary(M4_food_germany_ITT_Robust))[3], coef(summary(M4_food_germany_ITT_Robust))[4], coef(summary(M4_food_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_food_germany_ITT_Robust ))[2,2], coef(summary(M4_food_germany_ITT_Robust))[3,2], coef(summary(M4_food_germany_ITT_Robust))[4,2], coef(summary(M4_food_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Meat/Fish Consumption"), 4)),
  country = "Germany")

mydf_M4_food_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M4_food_china_ITT_Robust))[2],coef(summary(M4_food_china_ITT_Robust))[3], coef(summary(M4_food_china_ITT_Robust))[4], coef(summary(M4_food_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_food_china_ITT_Robust ))[2,2], coef(summary(M4_food_china_ITT_Robust))[3,2], coef(summary(M4_food_china_ITT_Robust))[4,2], coef(summary(M4_food_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Meat/Fish Consumption"), 4)),
  country = "China")


b4_Robust <- rbind(mydf_M4_food_china_ITT_Robust,mydf_M4_food_germany_ITT_Robust, mydf_M4_food_usa_ITT_Robust)

b4_Robust <- b4_Robust[, names(b4_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 


#Mobility

M4_mobility_usa_ITT_Robust <- lm_robust(slider_reduce_per_caruse_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="USA"  ))
summary(M4_mobility_usa_ITT_Robust)

M4_mobility_germany_ITT_Robust <- lm_robust(slider_reduce_per_caruse_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="Germany" ))
summary(M4_mobility_germany_ITT_Robust)

M4_mobility_china_ITT_Robust <- lm_robust(slider_reduce_per_caruse_z ~ treatment_mobility , data= subset (df, group =="Mobility"& country=="China" ))
summary(M4_mobility_china_ITT_Robust)

#Preparing Plots for Main Effects

mydf_M4_mobility_usa_ITT_Robust <- data.frame(
  coef = c(coef(summary(M4_mobility_usa_ITT_Robust))[2],coef(summary(M4_mobility_usa_ITT_Robust))[3], coef(summary(M4_mobility_usa_ITT_Robust))[4], coef(summary(M4_mobility_usa_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_mobility_usa_ITT_Robust ))[2,2], coef(summary(M4_mobility_usa_ITT_Robust))[3,2], coef(summary(M4_mobility_usa_ITT_Robust))[4,2], coef(summary(M4_mobility_usa_ITT_Robust))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Car Use"), 4)),
  country = "USA")

mydf_M4_mobility_germany_ITT_Robust <- data.frame(
  coef = c(coef(summary(M4_mobility_germany_ITT_Robust))[2],coef(summary(M4_mobility_germany_ITT_Robust))[3], coef(summary(M4_mobility_germany_ITT_Robust))[4], coef(summary(M4_mobility_germany_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_mobility_germany_ITT_Robust ))[2,2], coef(summary(M4_mobility_germany_ITT_Robust))[3,2], coef(summary(M4_mobility_germany_ITT_Robust))[4,2], coef(summary(M4_mobility_germany_ITT_Robust))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Car Use"), 4)),
  country = "Germany")

mydf_M4_mobility_china_ITT_Robust <- data.frame(
  coef = c(coef(summary(M4_mobility_china_ITT_Robust))[2],coef(summary(M4_mobility_china_ITT_Robust))[3], coef(summary(M4_mobility_china_ITT_Robust))[4], coef(summary(M4_mobility_china_ITT_Robust ))[5]),
  Treatment = factor(c("Animal\nWelfare", "Global\nClimate", "Local\nEnvironment", "Personal\nHealth")),
  se = c(coef(summary(M4_mobility_china_ITT_Robust ))[2,2], coef(summary(M4_mobility_china_ITT_Robust))[3,2], coef(summary(M4_mobility_china_ITT_Robust))[4,2], coef(summary(M4_mobility_china_ITT_Robust))[5,2]),
  DV = factor(rep(c("Intentions to Reduce Car Use"), 4)),
  country = "China")


c4_Robust <- rbind(mydf_M4_mobility_china_ITT_Robust,mydf_M4_mobility_germany_ITT_Robust, mydf_M4_mobility_usa_ITT_Robust)
c4_Robust<- c4_Robust[, names(c4_Robust) %in% c('Treatment', 'country', 'coef', 'DV', 'se')] 

d4_Robust <- rbind (c4_Robust,b4_Robust)




################ Final Figures ######################

#Figure 2

d1$type <- "OLS"
d1$ymin <- d1$coef - (1.96*d1$se)
d1$ymax <- d1$coef + (1.96*d1$se)

d1sub <- d1 %>% select(Treatment, coef, ymin, ymax, DV, country, type)

sum1_df$type <- "LASSOplus"
sum1_df$coef <- sum1_df$`50%`
sum1_df$ymin <- sum1_df$`2.5%`
sum1_df$ymax <- sum1_df$`97.5%`

sum1_dfsub <- sum1_df %>% select(Treatment, coef, ymin, ymax, DV, country, type)

e1 <- rbind(d1sub,sum1_dfsub)


Fig2 <- ggplot(e1, aes(x = Treatment, color = type, group =type, y =coef)) + 
  geom_errorbar(aes(ymin=ymin, ymax=ymax), size =2, width=.2, position=position_dodge(width =0.6)) +
  #geom_line(position=pd, size =2) +
  geom_point(aes (shape = type), position=position_dodge(width =0.6), size=6, fill="white")+ # 21 is filled circle +
  geom_hline(yintercept = 0, lty="dashed", size =2)+
  facet_grid(country~DV, scales ="free")+
  #scale_y_continuous(breaks=c(-.2,0,.2,.4,.6))+
  ylim(c(-0.3,0.8))+
  xlab("") +
  ylab("Estimated Effects") +
  #ggtitle("Classical Linear Regression vs. LASSOplus:\nFraming Effects on Concern about Car Use and Meat/Fish Consumption") +
  theme_bw() +
  theme(legend.justification=c(1,0),
        legend.position="bottom")      + 
  theme(plot.title = element_text(size = 34, face = "bold"),legend.title =element_blank(), legend.text = element_text(size = 30, face = "bold"),axis.text.x =element_text(size=30, face ="bold"),axis.text.y =element_text(size=30, face ="bold"),strip.text.x = element_text(size=30, face ="bold"),strip.text.y = element_text(size=30, face ="bold"), axis.title.y = element_text(size = 30, face ="bold"))+
  scale_colour_manual("type", values = c("LASSOplus" = "darkblue", "OLS" = "darkred")
  )+guides(fill=FALSE)


ggsave("Fig2.png", Fig2, height = 12, width = 11*2, dpi = 300)

#Figure 3

d2$type <- "OLS"
d2$ymin <- d2$coef - 1.96*d2$se
d2$ymax <- d2$coef + (1.96*d2$se)

d2sub <- d2 %>% select(Treatment, coef, ymin, ymax, DV, country, type)


sum2_df$type <- "LASSOplus"
sum2_df$coef <- sum2_df$`50%`
sum2_df$ymin <- sum2_df$`2.5%`
sum2_df$ymax <- sum2_df$`97.5%`

sum2_dfsub <- sum2_df %>% select(Treatment, coef, ymin, ymax, DV, country, type)

e2 <- rbind(d2sub,sum2_dfsub)

Fig3 <- ggplot(e2, aes(x = Treatment, color = type, group =type, y =coef)) + 
  geom_errorbar(aes(ymin=ymin, ymax=ymax), size =2, width=.2, position=position_dodge(width =0.6)) +
  #geom_line(position=pd, size =2) +
  geom_point(aes (shape = type), position=position_dodge(width =0.6), size=6, fill="white")+ # 21 is filled circle +
  geom_hline(yintercept = 0, lty="dashed", size =2)+
  facet_grid(country~DV, scales ="free")+
  #scale_y_continuous(breaks=c(-.2,0,.2,.4,.6))+
  ylim(c(-0.3,0.8))+
  xlab("") +
  ylab("Estimated Effects") +
  #ggtitle("Classical Linear Regression vs. LASSOplus:\nFraming Effects on Support for Policies to Reduce Car Use and Meat/Fish Consumption") +
  theme_bw() +
  theme(legend.justification=c(1,0),
        legend.position="bottom")      + 
  theme(plot.title = element_text(size = 34, face = "bold"),legend.title =element_blank(), legend.text = element_text(size = 30, face = "bold"),axis.text.x =element_text(size=30, face ="bold"),axis.text.y =element_text(size=30, face ="bold"),strip.text.x = element_text(size=30, face ="bold"),strip.text.y = element_text(size=30, face ="bold"), axis.title.y = element_text(size = 30, face ="bold"))+
  scale_colour_manual("type", values = c("LASSOplus" = "darkblue", "OLS" = "darkred")
  )+guides(fill=FALSE)

ggsave("Fig3.png", Fig3, height = 12, width = 11*2, dpi = 300)

#Figure 4

d3$type <- "OLS"
d3$ymin <- d3$coef - (1.96*d3$se)
d3$ymax <- d3$coef + (1.96*d3$se)

d3sub <- d3 %>% select(Treatment, coef, ymin, ymax, DV, country, type)

sum3_df$type <- "LASSOplus"
sum3_df$coef <- sum3_df$`50%`
sum3_df$ymin <- sum3_df$`2.5%`
sum3_df$ymax <- sum3_df$`97.5%`

sum3_dfsub <- sum3_df %>% select(Treatment, coef, ymin, ymax, DV, country, type)

e3 <- rbind(d3sub,sum3_dfsub)


Fig4 <- ggplot(e3, aes(x = Treatment, color = type, group =type, y =coef)) + 
  geom_errorbar(aes(ymin=ymin, ymax=ymax), size =2, width=.2, position=position_dodge(width =0.6)) +
  #geom_line(position=pd, size =2) +
  geom_point(aes (shape = type), position=position_dodge(width =0.6), size=6, fill="white")+ # 21 is filled circle +
  geom_hline(yintercept = 0, lty="dashed", size =2)+
  facet_grid(country~DV, scales ="free")+
  #scale_y_continuous(breaks=c(-.2,0,.2,.4,.6))+
  ylim(c(-0.3,0.8))+
  xlab("") +
  ylab("Estimated Effects") +
  #ggtitle("Classical Linear Regression vs. LASSOplus:\nFraming Effects on Willingness to Pay more for Motor Fuels and Meat/Fish Consumption") +
  theme_bw() +
  theme(legend.justification=c(1,0),
        legend.position="bottom")      + 
  theme(plot.title = element_text(size = 34, face = "bold"),legend.title =element_blank(), legend.text = element_text(size = 30, face = "bold"),axis.text.x =element_text(size=30, face ="bold"),axis.text.y =element_text(size=30, face ="bold"),strip.text.x = element_text(size=30, face ="bold"),strip.text.y = element_text(size=30, face ="bold"), axis.title.y = element_text(size = 30, face ="bold"))+
  scale_colour_manual("type", values = c("LASSOplus" = "darkblue", "OLS" = "darkred")
  )+guides(fill=FALSE)

ggsave("Fig4.png", Fig4, height = 12, width = 11*2, dpi = 300)

#Figure 5

d4$type <- "OLS"
d4$ymin <- d4$coef - (1.96*d4$se)
d4$ymax <- d4$coef + (1.96*d4$se)

d4sub <- d4 %>% select(Treatment, coef, ymin, ymax, DV, country, type)

sum4_df$type <- "LASSOplus"
sum4_df$coef <- sum4_df$`50%`
sum4_df$ymin <- sum4_df$`2.5%`
sum4_df$ymax <- sum4_df$`97.5%`

sum4_dfsub <- sum4_df %>% select(Treatment, coef, ymin, ymax, DV, country, type)

e4 <- rbind(d4sub,sum4_dfsub)


Fig5 <- ggplot(e4, aes(x = Treatment, color = type, group =type, y =coef)) + 
  geom_errorbar(aes(ymin=ymin, ymax=ymax), size =2, width=.2, position=position_dodge(width =0.6)) +
  #geom_line(position=pd, size =2) +
  geom_point(aes (shape = type), position=position_dodge(width =0.6), size=6, fill="white")+ # 21 is filled circle +
  geom_hline(yintercept = 0, lty="dashed", size =2)+
  facet_grid(country~DV, scales ="free")+
  #scale_y_continuous(breaks=c(-.2,0,.2,.4,.6))+
  ylim(c(-0.3,0.8))+
  xlab("") +
  ylab("Estimated Effects") +
  #ggtitle("Classical Linear Regression vs. LASSOplus:\nFraming Effects on Intentions to Reduce Car Use and Meat Consumption") +
  theme_bw() +
  theme(legend.justification=c(1,0),
        legend.position="bottom")      + 
  theme(plot.title = element_text(size = 34, face = "bold"),legend.title =element_blank(), legend.text = element_text(size = 30, face = "bold"),axis.text.x =element_text(size=30, face ="bold"),axis.text.y =element_text(size=30, face ="bold"),strip.text.x = element_text(size=30, face ="bold"),strip.text.y = element_text(size=30, face ="bold"), axis.title.y = element_text(size = 30, face ="bold"))+
  scale_colour_manual("type", values = c("LASSOplus" = "darkblue", "OLS" = "darkred")
  )+guides(fill=FALSE)


ggsave("Fig5.png", Fig5, height = 12, width = 11*2, dpi = 300)



################ Regression Tables ######################
#install.packages("stargazer")
library(stargazer)

#### Main OLS Results

#################Concern OLS

tablefig2a <- texreg::htmlreg(list(M1_food_china_ITT, M1_food_germany_ITT, M1_food_usa_ITT,M1_mobility_china_ITT,M1_mobility_germany_ITT, M1_mobility_usa_ITT),file="tablefig2a.html",
                        single.row = TRUE,
                        caption = "Table corresponding to Figure 2a: Concern", 
                        stars = c(0.01, 0.05),
                        digits = 3,
                        include.ci = F,
                        omit.coef = "(Intercept)",
                        custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                        custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                        #object.names = FALSE,
                        #notes.append = FALSE, notes.align = "l",
                        custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")


#################Policy Support

tablefig3a <- texreg::htmlreg(list(M2_food_china_ITT, M2_food_germany_ITT, M2_food_usa_ITT,M2_mobility_china_ITT,M2_mobility_germany_ITT, M2_mobility_usa_ITT),file="tablefig3a.html",
                              single.row = TRUE,
                              caption = "Table corresponding to Figure 3a: Policy Support", 
                              stars = c(0.01, 0.05),
                              digits = 3,
                              include.ci = F,
                              omit.coef = "(Intercept)",
                              custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                              custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                              #object.names = FALSE,
                              #notes.append = FALSE, notes.align = "l",
                              custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")

                        
#################Willingness to Pay

tablefig4a <- texreg::htmlreg(list(M3_food_china_ITT, M3_food_germany_ITT, M3_food_usa_ITT,M3_mobility_china_ITT,M3_mobility_germany_ITT, M3_mobility_usa_ITT),file="tablefig4a.html",
                              single.row = TRUE,
                              caption = "Table corresponding to Figure 4a: Willingness to Pay", 
                              stars = c(0.01, 0.05),
                              digits = 3,
                              include.ci = F,
                              omit.coef = "(Intercept)",
                              custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                              custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                              #object.names = FALSE,
                              #notes.append = FALSE, notes.align = "l",
                              custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")

#################Behavioral Intentions

tablefig5a <- texreg::htmlreg(list(M4_food_china_ITT, M4_food_germany_ITT, M4_food_usa_ITT,M4_mobility_china_ITT,M4_mobility_germany_ITT, M4_mobility_usa_ITT),file="tablefig5a.html",
                              single.row = TRUE,
                              caption = "Table corresponding to Figure 5a: Behavioral Intentions", 
                              stars = c(0.01, 0.05),
                              digits = 3,
                              include.ci = F,
                              omit.coef = ("Intercept"),
                              custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                              custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                              #object.names = FALSE,
                              #notes.append = FALSE, notes.align = "l",
                              custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")

#### Robustness Check OLS Results

#################Concern 

tablefig2a_Robust <- texreg::htmlreg(list(M1_food_china_ITT_Robust, M1_food_germany_ITT_Robust, M1_food_usa_ITT_Robust,M1_mobility_china_ITT_Robust,M1_mobility_germany_ITT_Robust, M1_mobility_usa_ITT_Robust),file="tablefig2a_Robust.html",
                              single.row = TRUE,
                              caption = "Table corresponding to Figure 2 Robustness Check: Concern", 
                              stars = c(0.01, 0.05),
                              digits = 3,
                              include.ci = F,
                              omit.coef = "(Intercept)",
                              custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                              custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                              #object.names = FALSE,
                              #notes.append = FALSE, notes.align = "l",
                              custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")


#################Policy Support

tablefig3a_Robust <- texreg::htmlreg(list(M2_food_china_ITT_Robust, M2_food_germany_ITT_Robust, M2_food_usa_ITT_Robust,M2_mobility_china_ITT_Robust,M2_mobility_germany_ITT_Robust, M2_mobility_usa_ITT_Robust),file="tablefig3a_Robust.html",
                              single.row = TRUE,
                              caption = "Table corresponding to Figure 3 Robustness Check: Policy Support", 
                              stars = c(0.01, 0.05),
                              digits = 3,
                              include.ci = F,
                              omit.coef = "(Intercept)",
                              custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                              custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                              #object.names = FALSE,
                              #notes.append = FALSE, notes.align = "l",
                              custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")


#################Willingness to Pay

tablefig4a_Robust <- texreg::htmlreg(list(M3_food_china_ITT_Robust, M3_food_germany_ITT_Robust, M3_food_usa_ITT_Robust,M3_mobility_china_ITT_Robust,M3_mobility_germany_ITT_Robust, M3_mobility_usa_ITT_Robust),file="tablefig4a_Robust.html",
                              single.row = TRUE,
                              caption = "Table corresponding to Figure 4 Robustness Check: Willingness to Pay", 
                              stars = c(0.01, 0.05),
                              digits = 3,
                              include.ci = F,
                              omit.coef = "(Intercept)",
                              custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                              custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                              #object.names = FALSE,
                              #notes.append = FALSE, notes.align = "l",
                              custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")

#################Behavioral Intentions

tablefig5a_Robust <- texreg::htmlreg(list(M4_food_china_ITT_Robust, M4_food_germany_ITT_Robust, M4_food_usa_ITT_Robust,M4_mobility_china_ITT_Robust,M4_mobility_germany_ITT_Robust, M4_mobility_usa_ITT_Robust),file="tablefig5a_Robust.html",
                              single.row = TRUE,
                              caption = "Table corresponding to Figure 5 Robustness Check: Behavioral Intentions", 
                              stars = c(0.01, 0.05),
                              digits = 3,
                              include.ci = F,
                              omit.coef = ("Intercept"),
                              custom.model.names = c("China: Meat/Fish","Germany: Meat/Fish","USA: Meat/Fish", "China: Car","Germany: Car","USA: Car"), 
                              custom.coef.names = c("Animal Welfare (Meat/Fish)", "Global Climate (Meat/Fish)", "Local Environment (Meat/Fish)","Personal Health (Meat/Fish)","Animal Welfare (Car)", "Global Climate (Car)", "Local Environment (Car)","Personal Health (Car)"),
                              #object.names = FALSE,
                              #notes.append = FALSE, notes.align = "l",
                              custom.note = "Table entries are standardised regression coefficients from an OLS regression with estimated robust standard errors in parantheses. ***p < 0.001, **p < 0.01, *p < 0.05, +p < 0.1")
